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ESTIMATION IN MIXED FREQUENCY DISTRIBUTIONS 


SUMMARY 


In the analysis of atmospheric data, distributions which are asymmetrical 
or multimodal are often encountered. These distributions are usually composed 
of two or more distinct homogeneous subpopulations and are designated as 
"mixed" or "compound" frequency distributions. 

The problem of estimating parameters in a mixed distribution is of 
considerably greater complexity than that of estimation in a single unimodal 
distribution. Not only is it necessary to estimate the parameters of each com- 
ponent, but the proportionality factors, which express the proportion or per- 
centage of each component in the mixture, must also be estimated. In the case 
of a mixture consisting of K components each having two parameters, there 
are 3K-1 parameters to be estimated in the resultant distribution. 

Karl Pearson solved such a problem as early as 1894 for a compound 
normal population, using the method of moments. The problem was investigated 
later, first by Charlier, and then in a joint effort by Charlier and Wicksell, 
who greatly simplified the theory. 

The papers presented in this document are concerned with estimation 
in compound distributions and are the results of an investigation performed by 
the Institute of Statistics, University of Georgia, Athens, Georgia, as a part 
of NASA contract NAS 8-11175 with the Terrestrial Environment Branch, 
Aerospace Environment Division, Aero-Astrodynamics Laboratory, NASA - 
George C. Marshall Space Flight Center, Huntsville, Alabama. Dr. A. C. 

Cohen, Jr. , is the principal investigator and the NASA contract monitors are 
Mr. O. E. Smith and Mr. L. W. Falls. 

The methods and procedures presented are practical and are applicable 
to experiments in which the data exhibit qualities that require a mixed distribu- 
tion as the statistical model. 



ESTIMATION IN A MIXTURE OF A POISSON WITH A 
NEGATIVE BINOMIAL DISTRIBUTION 

by A. Clifford Cohen, Jr. 


Introduction 


Many of the distributions encountered in the analysis of atmospheric 
data are the result of mixing two or more separate component distributions. 
Such distributions are therefore of particular interest to aerospace scientists. 
Mixtures of two Poisson distributions , mixtures of two exponential distribu- 
tions and mixtures of two normal distributions have been considered in 
previous memoranda [ 1,2,3]. This paper is concerned with estimation of 
the parameters in a distribution consisting of a Poisson component mixed with 
a negative binomial component. 

The Probability Function 

The density function of a mixed distribution with components fj(x) and 
f 2 (x) combined in the proportions a and ( 1 - a) respectively may be written 
as 


f(x) = af^x) + (1 - a) f 2 (x) . 


( 1 ) 


Let fj(x) be the density function of the Poisson component, which we 
write as 


-A x 

fl(x) = ^“^7 — 5 x = 0, 1, 2 . . . 


( 2 ) 


Let f 2 (x) be the density function of the negative binomial component, which 
may be written as 


f(x) = 


T(x + k) 
x! T(k) 


k. , . x 

P (1 - P) 


x = 0, 1, 2 . . . , 


(3) 


3 


L 



as 


f(x) = 


r(x + k) 
x! r(k) 


x 


(1 +q) 


k+x 


; x = 0, 1, 2 ... , 


(4) 


or as 


f(x) = 


r(x + k) 
r(k) 


x 


H)"W' 


x = 0, 1, 2 


(5) 


where q = ( 1 - p) /p , m = qk; m , k, q>0;0<p<i. 


Estimation of Parameters 


When f 2 ( x) assumes the form given in equation (5) , the first four 
factorial moments and the zero probability of the mixed distribution are 


= OiX + (1 - Q!)m, 

^[2] = ^ + ( 1 - “) (k + i)m 2 /k » 

= aA 3 + (1 - a) (k + 1) (k + 2)m 3 /k 2 , ► 

/x j = aA 4 + (1 - o?)(k + l)(k + 2)(k+ 3)m 4 /k 3 , 
f ( 0) = P = ae ^ + ( 1 - a) ( i + m/k) ^ . 


( 6 ) 


On equating the above moments to corresponding sample moments and 
on setting P = n 0 /n for a random sample of size n where n 0 is the number of 
zero observations in the sample, we obtain a system of equations, any four of 
which might be solved simultaneously for estimates of the parameters o;, 

X, m and k. Here, we consider estimates based on (a) the first four moments, 
and (b) on the first three moments and selected frequencies. 
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Equate M [ir M [2] and M [3] to x, p [ 2] and i> [3] 
first three moments of ( 6) become 


respectively; the 


x = aX + (1 - Q')m, 


= aX 2 + ( 1 - a) (k + i)m 2 /k , 


=aX 3 +(l-a)(k+l)(k + 2)m7k 2 , ) 


( 7 ) 


where 

R 

x= V x a /n , 

U r, X 

x=0 

R 

• i ~ V x(x - 1) . . . (x - j + l)n /n , 

lJJ x=o x ; 


( 8 ) 


in which n denotes the number of sample observations for which the random 
x 

variable X = x, and R denotes the largest observed value of X. 

Let us assume that X is known. We then solve the first equation of (7) 
for m and thereby obtain 


m = 


x - aX 
i - a 


(9) 


On substituting this value into the last two equations of (7) , and simplifying, 
we obtain 
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I 

II III! II II I I I I I Mil I II I I I II III I I III II I II I I I II I I I I I I I I I I I I II I I I I 1 1 I II I I I II III I I ■ I I 



(i - O') 


( o; | X) , 



( "t21 ~ rf> 
(x - aX) z 


\ 



(i - O') - o;X 3 ) 

(x - aX) (^j 2 ] " aX ^ 


ip2 ( a l X) . 


> 


/ 


( 10 ) 


With X thus assumed to be known, we eliminate k between the two equations of 
(10) to obtain a quadratic equation in a which we write as 


O! 2 - Bo: + C = 0 , 


(11) 


where 


(x*\ ol - 2i> 2 r , ) + X{v.. + 2xv , , ) - A 2 ( 4r r , - 3?) + X^x \ 


B = 


f 3] [2] [3] [2] 


12 ) 


A "[3] - 3X2 "t2] + 3X35E - 


c = 


x ^ [3] + (3? 2 ^ [2] )^ [2 ] 

^[3] " 3A 2 ^ [2] +3A 3 x- A 4 


( 12 ) 


The required estimate a * must be a positive root of (11) such that 
0 < a* <1 and such that corresponding estimates m* and k* are both positive. 
For some combinations of sample data, such an estimate may not exist and for 
some samples, both roots of (11) may satisfy the criteria set forth above for 
acceptable estimates of a. In the former case we can conclude only that either the 
sample data fail to conform to the model specified by equation ( 1) or the sample 
is too small. When the solution of equation (11) produces two acceptable 
estimates of a, we might choose the one which produces the closest agreement 
between the expected and the observed fourth moment or between the expected 
and observed number of observations in some specified class. For example, 
the zero class or perhaps the modal class might be chosen. 
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With a* determined, we estimate m from (9) as 


m 


* 


x - Acn* 
1 - a* 


( 13 ) 


and k from the first equation of (10) as 


k* = | X) - l] -1 , 


(14) 


or as 


V J2] " a * X * 

(1 - a *)( m *) 2 


(15) 


The asterisks employed here serve to distinguish estimators or estimates 
from the parameters being estimated. When desired, estimates of p and q 
as employed in (3) and (4) may be expressed as 


q* =m*/k* , P* = 1/(1 +q*) . 


(16) 


Estimates Based on First Four Moments . The foregoing results are 
based on the assumption that X is known. If indeed this assumption were true, 
our task would now be complete. In the more general case, where X too must 
be estimated from the sample data, we might begin by assuming a value X^ 


as a first approximation and subsequently compute first approximations ^ , 

m^ and k^ using (11), (13) and (14), or (15). On substituting these 

values into the fourth equation of ( 6) , we calculate a first approximation to 
the fourth factorial moment, n* for comparison with the corresponding 

1 J (l) 

sample moment v. . Once we have found two values X... and X,. .. in a 
[4] (i) (l+l) 

sufficiently narrow interval such that v ^ is in the interval (Mj 4 j ^ > 

P[ 4 ] ^ ) , we can interpolate linearly as indicated in Table I to determine 

final estimates a* , X* , m* , k* for which the first four sample moments are 
in agreement with the first four distribution moments. 
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TABLE I. LINEAR INTERPOLATION FOR FINAL ESTIMATES 


X 

a 

m 

k 

M [4] 

A (i) 

a a) 

m (i) 

k «) 

M [4](i) 

X* 

a* 

% 

m 

k* 

V U] 

X (i+i) 

a. . .. 

(i+l) 

m (i+D 

Vd 

M [4] (i+l) 


Estimates Based on First Three Moments Plus Frequency in a Selected 
Class. Because of larger sampling errors that are inherent in higher moments, 
it might sometimes be desirable to exchange the fourth equation of ( 6) which 
involved the fourth moment for an equation setting the observed proportion of 
observations equal to the expected proportion in a selected class. For example, 
the zero class might be chosen in which case we would employ the fifth equation 
of (6) rather than the fourth. In general, however, we might wish to employ 
the modal class. We begin with a first approximation and subsequently 


determine , m ^ 


and ^ from the first three moments using equations 


( ii) , (13) and ( 14) or ( 15) , as before. We then compute f^ ^ (x) for the 

selected value of x using equation (1) and compare the value thus obtained 

Subsequently, we determine second 


with the observed proportion n^/n 
approximations A^ , “(2)’ m (2) 


and k. 


in ( 1) , we compute f , 0 . (x) 

\ 6 ) 


( 2 )* 

As soon as two values A 


With these values substituted 


a) a “ d Vd in a 


sufficiently narrow interval have been found such that n^/n is in the interval 
[f^ (x) , f^ + ^ (x) ] , we interpolate linearly as indicated in Table 13 for the 
final estimates X** , a 


**, m* * , and k**. 


The double asterisks distinguish estimates based on the first three 
moments and a selected density from those based on the first four moments 
(designated with single asterisks) and from parameters being estimated 
(without asterisks) . 
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TABLE II. LINEAR INTERPOLATIONS FOR FINAL ESTIMATES 


X 

a 

m 

k 

f(x) 

x (i) 

“(i) 

m a) 

No 

V x) 

X** 

Of 5 ®* ^ 

m** 

k* * 

n /n 

X 

Vi> 

"(i+l) 

m (i+l) 

k (i+i) 

'« + i) (x) 


In both of the estimating procedures presented here, a first approxima- 
tion to X in the vicinity of x should usually be satisfactory. 


An Illustrative Example 


To illustrate the practical application of estimating procedures developed 
here, we construct a sample consisting of 1000 observations from a mixed 
Poisson and negative binomial distribution for which a = 0. 4, X = 3, m = 4. 5 
and k = 3. Data for this sample are shown in Table in. 


For this sample, n = 1000, x = 3. 915, ^ 2 ] = 0^4, ^[ 3 ] = *^4 . 946 

and = 1, 138. 536. On substituting these values into ( 12) , we have 

-274. 408722 + 291. 81222 X - 64. 808775 X 2 + 3. 915 X 3 


B = 


134. 946 X - 60. 102 X 2 + 11. 745 X 3 - X 4 


32. 65690365 

134. 946 X - 60. 102 X 2 + 11. 745 X 3 - X 4 ’ 
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TABLE HI. SAMPLE FROM MIXED POISSON AND 
NEGATIVE BINOMIAL DISTRIBUTION 



substituted in the above expressions for B and C, equation ( li) becomes 


a 2 - 1. 28688692 a + 0. 33881630 = 0 . 


Solving with the aid of the quadratic formula, we find 


a = ~ II, 28688692 ± *4 0. 30081276 ] 
i 1) 2 


from which 


a ^ = 0. 36921 


or 


0. 91768 . 
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When the smaller of these values is substituted into ( 13) and ( 15) , we 

find 


“(I) 


3. 915 - ( 0 . 36921) (3. 2 ) 
0. 63079 


4. 3335 , 


k (i) = I 1 - 37 2 ° 78 94 2 - l]" 1 = 2. 6876 . 


When the larger root of the quadratic equation is substituted into ( 13) 
and (15) , we obtain a negative value for k. Since, however, k must be 
positive, the larger root is rejected and the values given above as based on 
the smaller root are accepted as first approximations to m and k and the 
smaller root itself is accepted as our first approximation to a. Using the 
values thus calculated, we employ the fourth equation of ( 6 ) to calculate as a 
first approximation to the fourth factorial moment 


M [4] = 1165. 32 , 


which is to be compared with the sample value ^ 4 j = 1138. 536. 

For a second approximation, we choose X. . = 2. 8 , and following the 

( 2 ) 

same procedure as with the first approximations, we find a = 0. 42645, 

( 2 ) 

m ( 2 ) = 4.7440, k^ =3.4127 and ^[ 4 ]^) = ^45. 72. We ultimately try a 

third approximation =2.7 for which = 1138.27. The results of 

these different approximations are summarized in Table IV and the tables 
following. 


In solving equation (11) for various values of X as shown above, the 
smaller root was chosen as the appropriate approximation to the estimate of 
a in each case as the larger root, even when otherwise acceptable, resulted 
in a greater disparity between observed and expected number of observations 
in the zero class. 
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TABLE IV. FIRST APPROXIMATIONS 


X 


2. 6 

a 2 

2.7 

a 2 

2. 8 

a 2 

3. 0 

a 2 

3.2 

a 2 


Estimating Equations 
1. 09214560 a + 0. 31012438 - 0 
1. 13280360 a + 0.31327564= 0 
1. 16985786a + 0. 31702728 = 0 
1. 23410203 a + 0. 32645477 = 0 
1. 28688692 a + 0. 33881630 = 0 


TABLE V. SECOND APPROXIMATION 


X 

a 

m 

k 

q 

P 

M [4] 

2. 6 

Imaginary 






2.7 

0. 47960 

5. 0347 

3. 9421 

1. 2772 

0. 4391 

1138. 27 

2. 8 

0. 42645 

4. 7440 

3. 4127 

1. 3901 

0. 4184 

1145. 72 

3. 0 

0. 38403 

4. 4855 

2. 9614 

1. 5146 

0. 3977 

1155. 96 

3.2 

0. 36922 

4. 3335 

2. 6876 

1. 6124 

0. 3828 

1165. 32 


Additional computations were made with X = 2. 6 and X = 3. With 
X = 2. 6, equation (11) has no real roots and thus this value must be rejected. 
Approximations corresponding to X = 3 are included in Table V. For final 
estimates based on the first four moments, we interpolate between approxima- 
tions corresponding to X = 2. 7 and X = 2. 8 as shown in Table VI. 
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TABLE VI. FINAL ESTIMATES BASED ON 
FIRST FOUR MOMENTS 


X 

a 

m 

k 

M [4] 

2 . 8 

0. 42645 

4. 7440 

3. 4127 

1145.72 

2. 704 

0. 47777 

5. 0242 

3. 9229 

1138. 54 

2. 7 

0 . 47960 

5. 0347 

3. 9421 

1138.27 


Accordingly, final moment estimates are X* = 2.704, a* = 0. 4778, 
m* = 5. 024, k* = 3. 923, q* = 1. 280, p* = 0. 439, which are to be compared 
with the actual distribution parameters, X = 3, a = 0. 4, m = 4. 5, k = 3, 
q = 1. 5, p = 0. 4. 

Somewhat improved estimates are obtained when the frequency in the one 
class rather than the fourth factorial moment is utilized as the basis for our fourth 
estimating equation. Estimates, thus based on the first three moments and the 
frequency in the one class are determined by interpolation as indicated in 
Table VII, where for a given value of X, f( i) is calculated using the probability 


TABLE VII. INTERPOLATION DETERMINING ESTIMATES BASED 
ON THE FIRST THREE MOMENTS 


A 

a 

m 

k 

f(D 

3. 0 

0. 38403 

4. 4855 

2. 9614 

0. 12896 

2. 996 

0. 38495 

4. 4911 

2. 9712 

0. 12900 

2 . 8 

0. 42645 

4. 7440 

3.4127 

0. 13080 


function ( 1 ) with values for a, m, and k based on the first three moments, 
determined as previously described. Final estimates based on the first three 
moments and the frequency of ones are thus X** = 2. 996, a** = 0. 385, 
m** = 4 . 49 ^ k** = 2 . 97 and from (16) q** = 1. 51 and p** =0. 40. 
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Essentially these same estimates are obtained, when our calculations are based 
on frequencies in the zero, two, three, or four classes or even on a combination 
of these classes rather than on the one class as has been done here. The im- 
provement achieved by equating observed frequencies to expected frequencies 
rather than equating the observed fourth moment to the expected fourth moment 
apparently reflects the lack of stability in the higher sample moments. The 
fact that very small changes in a result in quite large variations in the higher 
moments is also involved. For comparison all estimates calculated are listed 
in Table Vin alongside the true parameter values. 


TABLE VIII. ALL ESTIMATES ALONGSIDE 
TRUE PARAMETER VALUES 


Basis 

X 

a 

m 

k 

q 

P 

Parameter Values 

3 

0. 4 

4. 5 

3 

1. 5 

0.4 

First Four Moments 

2. 704 

0. 478 

5. 02 

3. 92 

1. 28 

0.44 

First Three Moments 

2. 996 

0. 385 

4. 49 

2. 97 

1. 51 

0. 40 

and Freq. of Ones 
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ESTIMATION IN MIXTURES OF TWO 
NORMAL DISTRIBUTIONS 

by A. Clifford Cohen, Jr. 

Summary 

This paper is concerned primarily with the method of moments in 
dissecting a mixture of two normal distributions. In the general case, with 
two means, two standard deviations, and a proportionality factor to be estimated, 
the first five sample moments are required, and it becomes necessary to find a 
particular solution of a ninth degree polynomial equation that was originally 
derived by Karl Pearson [ 4 ] . A procedure which circumvents solution of the 
nonic equation and thereby considerably reduces the total computational effort 
otherwise required, is presented. Estimates obtained in the simpler special 
case in which the two standard deviations are assumed to be equal are employed 
as first approximations in an iterative method for simultaneously solving the 
basic system of moment equations applicable in the more general case in which 
the two standard deviations are unequal. Conditional maximum likelihood and 
conditional minimum chi-square estimation subject to having the first four 
sample moments equated to corresponding population moments, are also con- 
sidered. An illustrative example is included. 

Research was sponsored by the Aero-Astrodynamics Laboratory of the 
Marshall Space Flight Center, National Aeronautics and Space Administration; 
Contract NAS8-11175. 


Introduction 


Distributions which result from the mixing of two or more component 
distributions are designated as "compound" or "mixed" distributions. They 
may be further described by designating distribution types of the individual 
components. Such distributions arise in a wide variety of practical situations 
ranging from distributions of wind velocities to distributions of physical 
dimensions of various mass produced items. Compound normal distributions 
were studied as early as 1894 by Karl Pearson [4] and later by Charlier [5] 
and by Charlier and Wicksell [6]. More recently, compound Poisson and 
compound exponential distributions have been studied by Rider [ 7 ] and by 
Cohen [1, 2,8,9] .Compound binomial distributions have been studied by 
Blischke [10]. 
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This paper is concerned with estimation of the parameters 6 j, crj, 0 2 , 
cr 2 , and a of the compound normal distribution with density 


f(x) = of i(x) + (1 - Q!)f 2 (x) , 


( 1 ) 


where 


fj(x) = (2 no\) ^ exp-[(x - 0 1 )/cr 1 ] 2 /2, 

> ( 2 ) 

f 2 (x) = (2ircr 2 ) -1 ^ exp-[(x - 0 2 )/o 2 ] 2 / 2 . 

j 

Karl Pearson [4] derived estimators for the parameters of this distribution by 
equating sample moments to corresponding population (theoretical) moments. 
The evaluation of his estimators involved the solution of a ninth degree poly- 
nomial equation. Before the advent of modern electronic computers, this was 
considered a rather formidable obstacle to the application of his results. 
Charlier and Wicksell [6] succeeded in considerably simplifying Pearson’s 
results. The present effort represents an attempt to simplify these estimators 
further in order that they might be more readily available for use in appropriate 
practical applications. A procedure is presented whereby the direct solution 
of Pearson's nonic equation can be circumvented through an iterative process 
which involves solving a cubic equation for a unique negative root. In addition 
to considering the most general case in which all five parameters of ( 1) 
must be estimated from the data available in a given sample, some of the 
more important special cases in which one or more of the parameters are 
known in advance of making sample observations are also examined. 


Estimation in the General Case 

Except for a few changes in notation, the presentation of this section 
is essentially that of Charlier and Wicksell [6]. 

The kth moment of f(x) taken about the origin may be written as 


u[ = a 


oo oo 

J x^fj(x)dx + ( 1 - a) J x^f 2 (x)dx , 


(3) 
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where fj(x) and f 2 (x) are density functions of the two component distributions 
as given in equations ( 2) . The kth central moment becomes 


M f (x - /xi) k f 1 (x)dx + (1 - a) J (x - pi) k f 2 (x)dx 


(4) 


If we let 


mj = - pj and m 2 = 0 2 - p| , 


(5) 


the first non-central and the first five central moments of (1) follow as 


Pi = <20! + (1 - a)e 2 , 


Pj = ami + ( 1 - a ) m 2 = 0 , 


P 2 = <*(oi + m^) + (1 - a) (a| + m^) , 


p 3 = am^ 2 ! + m 2 i) + (1 - a)m 2 (3cr 2 + m|) , 


p 4 = a [ 3cr * + 6m 2 cr 2 + mf ] + ( 1 - a) [ 3a 3 + 6m 2 o| + m 2 ] , 

p 5 = amj [I5cr j + 10m 2 cr \ + m j] + ( 1 - a) m 2 [ l5a| + 10m| o| + m 2 ] . 


> ( 6 ) 


Without any loss of generality let us suppose that ^ 0 2 . Accordingly 
< pj < 0 2 > and m i — 0 < m 2 . 


Upon equating population moments to corresponding sample moments , 
it follows from ( 6) that 


a[0 1 - x] + (1 - a) [0 2 - x] = 0 , 


(7) 
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and further that 


> 

ami + (1 - a)m 2 = 0, 

a[Oi + mi - r 2 ] + (1 - a) [o| + m 2 - ^ 2 ] = 0, 

a [3o^mi + m| - v 3 ] + ( 1 - a) [3ofm 2 + m 2 - r 3 ] =0, 

a [ 3af + 6miOi + mf - r 4 ] ' ( 8) 

+ ( 1 - a) [ 3a 2 + 6m 2 CT| + m 2 - i> 4 ] = 0 , 


a [ I5crimi + lOOimi + mf - ^ 5 ] 

+ ( 1 - a) [l5o 2 m 2 + lOofmf + mf - t'g] = 0, 


where x is the sample mean and (i = 2, 3, 4, . . . ) is the ith central 
moment of the sample. Equations (8) accordingly constitute a system of five 
equations to be solved simultaneously for estimates of the five parameters 
a, m 1} Oi, m 2 , a 2 . 

We eliminate a between the first and subsequent equations of ( 8) in 
turn and thereby reduce this system to the following four equations in the four 
unknowns <Ti, mi, ff 2 , m 2 : 


(Tj + mj - r 2 _ q-| + m| - r 2 


m i 


m 2 


P , 


3o jmj + mj - r 3 _ 3cr 2 m 2 + m 2 - r 3 
mi m 2 


\ 


( 9 ) 


/ 


( Equation ( 9) concluded on following page. ) 
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3 cr\ + 6m + m j - _ 3crf + 6m 2 (7 2 + m 2 - Vj 


2„2 


mi 


m 2 


l5cr^m 1 + lOa^mi + mj - y 5 _ l5cr|m 2 + 10g| m 2 + m a ~ 


mi 


m 2 


(9) 

(Con- 

cluded) 


In view of the presence of m A and m 2 in the denominators of (9) , we recognize 
that these equations and those subsequently derived from them are not valid in 
the symmetric case to be dealt with later in which m A = m 2 = 0. In all subse- 
quent results presented in this section it is understood that mi < 0 and 
m 2 > 0. 


With the introduction of in the first equation of ( 9) it follows that 
= m^ + v 2 - m\ , I 

i ( 10 ) 

cr 2 = m 2 /3 + v 2 - mj . I 

On replacing cr\ and <r 2 in the second, third and fourth equations of (9) 
with the values given in ( 10) above, it follows after considerable algebraic 
manipulation that 


m 1 m 2 [3/3 - 2(mj + m 2 )] = - v 3 , 

3/3 2 - 2(m 2 + itijirij + m 2 ) ] = -k 4 , 

► 

mim 2 [ 15(mi + m 2 )(3 2 - 20(mi + mjm 2 + m 2 )/3 
+ 6(m! + m 2 ) (m 2 i + m|) ] = -k 5 , 

/ 


(ID 
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where k 4 and k 5 are respectively the fourth and fifth order sample cumulants 
or semi-invariants; i. e. , 


k 4 = r 4 - Zv\ , ) 


( 12 ) 


k 5 = r 5 - iOv 2 v 3 • J 

When referring to population (theoretical) cumulants, we employ the Greek 
kappa thus: 


k 4 ~ A*4 " > K 5~ ^5 ~ 10^2^3 • 


Equations (11) accordingly constitute a system of three equations in the 
three unknowns, j3, mj and m 2 . 


In order to further simplify the system of equations (11) , let 


r = m j + m 2 , and v = m 4 m 2 


(13) 


When these transformations are introduced into equations (11), the system 
becomes 


3/3v - 2rv = -v 3 , ^ 


3/3 2 v - 2v(r 2 - v) = -k 4 , 


> 


(14) 


l5vr/3 2 - 20v(r 2 - v)/3 + 6vr(r 2 - 2v) = -k 5 . / 


On making the further transformation 


w = rv , 


(15) 


20 



the system of equations ( 14) becomes 


3/3v - 2w = -v 3 , 

3/3 2 v 2 - 2w 2 + 2v 3 = -k 4 v , 

1 5v 2 w/3 2 - 20vw 2 /3 + 20v 4 /3 + 6w 3 - 12v 3 w = -k 5 v 2 . , 


( 16 ) 


We now eliminate /? between the first and the second and between the 
first and the third equations of (16) , and our system becomes 

2w 2 - 6 v 3 - 3vk 4 + 4wi^ 3 - v\ = 0 , 

5w(2w - ^ 3 ) 2 - 20w 2 (2w - v 3 ) / (17) 

+ 20v 3 (2w - ^3) + 18w 3 - 36wv 3 = -3k 5 v 2 . I 


By introducing the further transformation 
z = w + ^3 , 

the two equations of ( 17) become 


2z 2 = 6v 3 + 3k 4 v + 3^1 , 


2z 2 ( z -r 3 ^ 3 ) + z(ul - 4v 3 ) + 3^1 + 24 v 3 v 3 = 3k 5 v 2 . j 


(19) 


When the expression for 2z 2 from the first of the above equations is substituted 
into the second of those equations , we obtain 
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( 20 ) 


-6^ 3 v 3 + 3k 5 v 2 + 9^ 3 k 4 + 6^ 3 
2v 3 + 3k 4 v + 4^1 


and when this value is reinserted into the first equation of ( 19) , we obtain a 
polynomial equation of the ninth degree in v, which for convenience we write 
as follows: 


a 9 v 9 + a 8 v 8 + ayv 7 + a 6 v 6 + a 5 v 5 + a 4 v 4 + a 3 v 3 
+ ajv 2 + a 4 v + a 0 = 0 , 


where 

a 9 - 8 , ' 

a 8 = 0 , 

a T = 28 k 4 , 

a 6 = 12i/| , 

a 5 = 30k 4 + 24k 5 i' 3 , 

a 4 = 148k 4 f 3 - 6k| , > 

a 3 = 96 v 3 - 36 y 3 k 4 k 5 + 9 k| , 

32 — -(21i^ 3 k 4 + 24^|kg) , 
a l = -32^fk 4 , 
a 0 = - 8v| , 


( 21 ) 


( 22 ) 


This is the well known nonic which was first given in 1894 by Pearson [4]. 

Since her© m 4 < 0 and m 2 > 0, then v = m 1 m 2 < 0, and the required - 
estimate v* is to be found as a negative real root of (21) . Throughout this 
paper asterisks (*) are employed to distinguish estimates from the parame- 
ters being estimated. Prior to the advent of electronic computers, the task 
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of solving this equation would have been considered a formidable assignment. 
Today, however, modern computers are available to perform the otherwise 
burdensome calculations involved. The ninth degree polynomial can be 
evaluated for any desired value of v in the vicinity of v* either by straight- 
forward substitution or by synthetic division. Standard iterative procedures 
will quickly lead to the required value of v* . Once v* is determined with 
the desired degree of accuracy, the estimate w* follows from (20) and (18) 
as 


-8»a,v* 3 + 3k K v* 2 + 6y a k A v* + 2v% 
2v* 3 + 3k4V* + 4t>| 


From ( 15) , we have r* = w*/ v* , and accordingly 


* - 8i> 3 v* 3 + 3knv* 2 + 6^ 3 k 4 v* + 

r " = - ■ v*T ar« + 3k 4 v* + ¥f) a 


(23) 


It follows from the defining relations ( 13) , that estimates of m 4 and 
m 2 are the roots of the quadratic equation 


Y 2 - r* Y + v* = 0 . 


(24) 


Thus, we have 

•< ' 

m* = t- [r* - \l r* 2 - 4v* }, 

A 

> 

m* = | [r* + \Ir* 2 - 4v* ] . 

Using the first equation of ( 14) , we estimate /3 as 


(25) 


j8* = | [2r* - v 3 /v* ] , 


(26) 
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and from (10) , we have 


of 2 = mf/3* + v 2 - mf , 

► 

of 2 = mjf/3* + - mf 2 . 

/ 

Finally, from (5) and from the second equation of ( 6) , we have 


0? = m* + x , 

02 * = m* + x , ► 

a* = - m* ) . ; 


(27) 


(28) 


Attention is again drawn to the fact that the preceding results are valid 
only if 0i# 0 2 . The symmetric case in which 0 4 = 0 2 and thus mj = m 2 = r = 
v = w = 0 is treated separately in the section entitled "Estimation in the 
Symmetric Cases. " 

Unfortunately, for some combinations of sample data the nonic 
equation (21) may have more than one negative root and accordingly we must 
choose between two or more sets of estimates. This lack of uniqueness bothered 
Pearson, and he suggested choosing the set of estimates which resulted in the 
closest agreement between the sixth central moment of the sample and the 
corresponding moment of the "fitted" compound curve. 

On Circumventing the Nonic Estimating Equation . In the event that 
r is known, we need only consider the first two equations of (14) , and when 
/3 is eliminated between them , we have the following cubic equation in v : 


6v 3 - 2r 2 v 2 + (3k 4 - 4r^ 3 )v + v\ = 0 . 


(29) 


Using Descartes' rule of signs, we find that this equation has one negative root 
plus either two positive or a pair of imaginary roots. The negative root is 
the required estimate v* . 
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Using this value for v* and the known value for r, the required esti- 
mates follow from (25), (26), (27) and (28) as before. 


Even though r is not known a priori, we might assume a value and 
employ the foregoing results to determine approximations to the required esti- 
mates which in turn can be substituted into the final equation of (6) , to approxi- 
mate the fifth central moment, 


Let r^j designate the ith approximation to r* and let designate 

the itii approximation to /if based on r^ . It should be relatively easy to find 

approximations r^ and r^ + ^ such that the sample moment v 5 is in the 

interval (u ... » M r /. As shown in the following section, a satisfactory 

5(i) 5(1+1) 

initial approximation to r can usually be found by assuming oj = <r 2 and solving 
the appropriate estimating equation for this special case. Once the interval 
between r^ and r^ + ^ has been narrowed sufficiently, the required esti- 
mate r* can be obtained by simple linear interpolation as indicated below. 


Ms 


M 


5(i) 


(i) 


r* v^_ 

r (i+l) M 5(i+1) 


With r specified, the well known method of Horner which utilizes 
synthetic division procedures is quite effective in solving (29) for v. Any 
standard iterative method, however, might be employed. A "trial and error" 
procedure based on linear interpolation with direct substitution in (29) , though 
perhaps not very economical of computational effort, will generally give 
satisfactory results. 

Various special cases in which one or more of the parameters are 
known or in some way restricted are sometimes of interest. With fewer 
parameters to be estimated, the number of sample moments involved is cor- 
respondingly reduced and the estimating equations are accordingly simpler 
and easier to solve. Some of the more important special cases are considered 
in the next two sections. 
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Conditional Maximum Like lihood and_ Conditional^ Minimum Chi-Square 
Estimation . In order to eliminate the effect of sampling errors resulting from 
direct use of the fifth order moment, let us consider a conditional maximum 
likelihood procedure. The first four sample moments are equated to corres- 
ponding population moments, and subject to this condition, r is determined so 
as to maximize 


n 

L r (r) = n f(x. |x, v 2 , v z , p 4 ) . 
i=l 1 


Since derivatives of L'(r) are somewhat unwieldy, the value of r which maxi- 
mizes L'(r) can conveniently be determined in most practical applications by 
actually calculating L'(r) for several values of r in the vicinity of its maxi- 
mum and employing either finite difference or graphical techniques. 

Grouped Data . In the case of grouped data, the conditional likelihood 
function might be expressed as 


L'(r) = II p 4 
i=l 


n i 


where k is the number of groups or classes into which the sample has been 

divided, n. is the number of sample observations in the ith class, x. is the 
1 — 1 

upper boundary of the ith class, and 


x. 




f(x|x, 


i-1 


V Z , v 4 )dx . 


In practice it is sometimes more convenient to minimize chi-square than to 
maximize L’(r). Kendall [11 (Vol. n, pp. 55-56)] has shown that with 
grouped data, the two methods are equivalent to the order n -1 /^ . We are 
therefore free to choose the method requiring the least computational effort. 
In the method of minimum chi-square, we seek the value of r which results 
in the minimum value for 
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k (n. - e.) 2 

X 2 (r|x, v 2 , v z , v A ) = 2 X e 1 

i=l i 


where ej = np^ is the expected number of observations in the ith class subject 
to the restriction that the first four sample moments are equated to corres- 
ponding population moments. As previously indicated, n is the number of 
sample observations in the ith class. 

In practice, chi-square can be calculated for several values of r in 
the vicinity of its minimum and the desired value of r can be determined 
either graphically or by employing finite difference techniques. 

Unrestricted maximum likelihood estimates would, of course, be 
preferable to any moment estimates, but with five parameters to be esti- 
mated, the unrestricted maximum likelihood estimating equations become 
quite intractable. The conditional maximum likelihood procedures suggested 
here are believed to represent a feasible compromise between the need for 
estimating equations that are easy to solve and the need for reliable estima- 
ting procedures. 


Estimation When a t = <r t = <r 

Here, we need only estimate the four parameters, 0*, 0 2 , ot and cr 
where as in the general case 0 A < 0 2 . Accordingly only the first four moments 
and or cumulants are required. Charlier and Wicksell [6] and Rao [12] 
among others have previously considered this special case. With 
0*1 = 0*2 = o*, equations (10) of the general five-parameter case now become 


a 2, = mjjS + v 2 - mi = m 2 /3 + ^ ~ m ! * 


(30) 


From this it follows that 


/ 3 =m 1 + m 2 = r , 

► 

o 2 - m 1 m 2 + v 2 = v + ^2 9 

* 

where r and v are as defined in ( 13) . 


(31) 
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The first two equations of ( 14) which are applicable here now become 


3rv - 2rv = -v 3 , 


3r 2 v - 2v(r 2 - v) = -k 4 , 


and subsequently 


rv = -v 3 , 

* 

T^v + 2v 2 = -k 4 . 

i 

From the first of the above equations 


(32) 


r = -v 3 /v . 


(33) 


When this value is substituted into the second equation of (32) , and both sides 
are multiplied by v, the applicable estimating equation becomes 


2v 3 + k 4 v + yjj = 0 . 


(34) 


From Descartes’ rule of signs it follows that unless v\ = 0, equation (34) has 
a single negative root, which is the required estimate v* regardless of whether 
k 4 is positive or negative. The other two roots are of no interest to us here. 

It is relatively easy to solve (34) for v* using standard iterative procedures, 
and from ( 33) it follows that 


r* = -^s/v* . 


(35) 
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With r* and v* thus determined, the estimates m* , m*, 0*, df 
and a* are given by (25) and (28) as in the general case, while ct * 2 follows 
from the second equation of (31) as 


cr* 2 =v* + 




(36) 


A first approximation to r i n the general case . In view of the relative 
ease with which estimates can be calculated when cr* = <t 2 , we are thus provided 
with a simple procedure for obtaining initial approximations to r in the general 
case. With <7j assumed equal to cr 2 , we obtain v^ as the negative root of 

(34) , and calculate r. . = -vJv. . , from (33). Unless the disparity between 

(o) (o) 

ctj and cr 2 is quite large, the resulting value, r^ , Provides a satisfactory 

starting point from which we can iterate to the final estimate r* in the general 
case. 


Estimation in the Symmetric Cases 

The compound normal distribution is symmetrical if (a) a = 1/2 with 
dq = cr 2 = <7, and if (b) 0 t = 0 2 . In the former instance the compound distribu- 
tion has been shown to be bimodal when 0 2 - 0j > 2cr, but otherwise unimodal 
[8]. In the second case, the resultant distribution is always unimodal. In the 
limiting (trivial) case in which 0j = 0 2 = p and cr 1 = cr 2 = cr, the resultant dis- 
tribution degenerates into a single normal distribution (p, cr) . 

Symmetric case with a = 1/2 and o~i = cr?, (61# 0?) . In this case, we 
need only estimate the three parameters cr, 6 t and 0 2 . From the second 
equation of (6) , with a = 1/2, it follows that mj = -m 2 and from the first 
equation of (31) , r = (3 = 0. Consequently, from the second equation of (32) , 
we have 


2v 2 + k 4 = 0 . 


(37) 


With the vanishing of the odd central population moments, equation (37) might 
have been obtained as a special case of (34) with v 3 replaced by zero, which 
in this instance is the appropriate estimate of p 3 . 
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It follows from ( 37) that 


v* = - \/-k 4 /2 . 


( 38 ) 


With v* given by (38) and with r = 0, the required estimates follow from 
(25) , (28) , and (36) as 


m* = n/-v* , mf = -m* , ^ 

0* = m* + x , 0* = mf + x , > 

cr* 2 = v * + v 2 . i 


(39) 


Symmetric case with 0i = 0? = 0 . Since estimating equations ( 9) in- 
volve division by m 4 and m 2 , and since it follows that m 4 = m 2 = 0 when 
0j = 0 2 , neither equations (9) nor subsequent equations derived from them 
are applicable here. We estimate 0 from (7) as 


0* = x . 


(40) 


With the vanishing of the odd central moments, estimation of the three remain- 
ing parameters , a, cr i and cr 2 necessitates use of the second, fourth, and 
sixth central moments. Applicable estimating equations accordingly become 


a(o\ - v 2 ) + (1 - a) {<4 - v 2 ) 


= 0 , 


\ 


C*(3cr \ - U A ) + ( i - a) (3cr| - v 4 ) =0, ? 


(41) 


a(15crf-i' e ) + ( 1 - a) ( 15af - v 6 ) = 0 . ) 
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The first two of these equations follow directly from (8) , and the third follows 
from (4) with k = 6 when fx e is estimated as v e . 

From the first equation of (41) , we have 


a = 


Ez. 

A 


- ° 2 
-A 


and (1 - a) = 


A 


- V 0 


A - A 


(42) 


On substituting these values for a and ( 1 - a) into the second and third equa- 
tions of ( 41) , after considerable algebraic manipulations we obtain 


> 

(A - v 2 ) (ct! - v 2 ) = -k4/3 , 

► 

(A - (°i - + A - = -kg/i5 , 


(43) 


where k 4 and k 6 are the fourth and sixth order sample eumulants respectively; 
i. e. , 


> 

k 4 = v 4 - 3^1 , 

► 

k 6 = v 6 - 15^2 - 10^1 + 30^1 . 

/ 

Let 


ti = cr 2 ! - v 2 and t 2 = cr| - v 2 , 


(44) 


(45) 


and equations ( 43) become 


t 4 t2 — kj3 , 


tjta ( t^ + t2) — kg/15 . 


\ 


> 


/ 


(46) 
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It follows from ( 46) that t* and are roots of the quadratic equation 

Y 2 - (ke/5k 4 )Y - (V3) = 0 . (47) 

Accordingly 

ti = | [ ( kg/ 5k 4 ) - \l ( k 6 / 5k 4 ) 1 + (4^3) ] , 

if = | I (ke/5k 4 ) + si ( kg/ 5k 4 ) 2 + (4V3) ] . 

With t* and tf thus determined, it follows from (45) and from (41) that 
of = t* + v 2 , of 2 = if + ^2 , a* = -tf/(t? - tf ) , (49) 

and of course 0* = x as given in ( 40) . 

Determining Which Case Is Applicable 

In the absence of a priori information concerning whether or not one of 
the symmetric special cases is applicable in lieu of the general case, the 
following criteria provide a basis for resolving this issue. 

(a) If = 0 and if k 4 < 0, the compound distribution is symmetric with 
CTj = cr 2 and with a = 1/2. 

(b) If = 0 and if k > 0, the compound distribution is symmetric with 

01 = o 2 . 

(c) If = 0 and /c 4 = 0, then 0 4 = 0 2 , cr i = <r 2 , and the "resultant!' 
distribution is in fact a single normal distribution. 

Of course the third central moment is zero in all symmetrical distri- 
butions, and the converse likewise is true. Therefore, the foregoing statements 
can be verified by examining applicable expressions for the fourth cumulant, 
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which is defined as K 4 = M 4 - 3/-c|- Using expressions from equations ( 6 ) for 
/U 4 and a< 2 » the fourth population cumulant in the most general case follows as 


k 4 = 3a(i - a) [(a 2 - <r 2 ) + (m 2 - m 2 ) ] 2 - 2[amf + (1 - ajmll . (50) 


When a = 1/2 and Qj = this implies that -m* = m 2 . With these values 
substituted into ( 50) , we have 


k 4 : -[mj + m 2 ] < 0 . ( 51) 

When 0! = 0 2 , this implies that m t = m 2 = 0, and in this instance (50) 
becomes 


k 4 = 3«(1 - a ) ( a , - ct 2 ) 2 > 0 . 


(52) 


When 6 < = 0 9 and a 1 = 0 - 9 , it follows from ( 52) that in this limiting 
case /c 4 = 0 , since (52) is applicable in all cases where 0 1 = 0 2 . 

In practical applications our classification problem is reduced to that 
of utilizing the sample statistics v 3 and k 4 = - 3v% in choosing the most 

acceptable hypothesis from among the following alternatives: H , 

Ulj H3 — U j U 

H i:ju 3 =0, k 4 <0’ H 2:p 3 =0,/c 4 >0’ H 3:M 3 *0' 

An Illustrative Example 

To illustrate the practical application of computational procedures de- 
veloped in this paper, we consider a mixed sample obtained by combining two 
separate component samples consisting of 334 and 672 observations respectively 
from normal populations. For these two component samples, x^ = 47. 716, 
s 2 ! = 33. 5663, N t = 334, x 2 = 57. 607, s| = 9. 1790 and N 2 = 672. For the 
resultant mixed sample, n = N* + N 2 = 1006, x = 54. 32306, v 3 = s 2 = 38. 9753, 
v 3 = -233. 876, = 5365. 13, v 5 = -67,821 and k 4 = 807. 91. In an effort to 

compensate for errors caused by grouping, Sheppard's corrections [H 
(Vol. I, p. 71)] have been applied both in computing s 2 and s 2 for the 
separate components and in computing moments of the resultant mixed sample. 
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We obtain a first approximation to v by substituting the values given 
above for k 4 and v z into equation (52) and solving for the negative root. By 
thus employing equation ( 52) , we implicitly assume = a 2 for this first 
approximation. In subsequent approximations , we abandon this restriction and 
accordingly replace equation (52) with (29) . For the first approximation our 
estimating equation, after making proper substitutions and simplifying, 
becomes 

V 3 + 403. 955 v + 27, 349 = 0 . 


With the aid of a desk calculator, straightforward substitution quickly 
yields as a solution the first approximation v^ ^ = -25. 7. A first approxima- 
tion to r then follows from ( 53) as 


r (i) = _( " 233 ' 876) /( -25. 7) = -9. 10 . 


Using these values, first approximations to the basic parameters follow from 
(54) , (25) and (28) as cr 2 (1) = = 14. 61, ™ 1(1) = -H- 36, ™ 2(1) = 2 - 2< 

and ot =0. 166. On substituting these values into the last equation of (6) , 


( 1 ) 

we subsequently calculate = -62,396 which is to be compared with the 

sample value, = -67, 821. 

As a second approximation to r, we let r^ = -4. 50 since further 

examination of the sample data with due regard for the shape of the histogram 
indicates that a value in the vicinity of -4 or -5 should be a good choice. This 
time we employ equation (29) rather than (34) in determining our new ap- 
proximation to v. With r = r , . = -4. 50, and with k 4 and v z as previously 

given, equation (29) becomes 


v 3 - 6. 75 v 2 - 297. 7 v + 9, 116 = 0 . 


On solving for the negative root of this equation, we find as our second 
approximation v^ = -23. 13. Using the above values for r^ and v^ , 
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equations (25), (26), (27), and (28) yield as second approximations to the 
remaining parameters of interest, m 


1 ( 2 ) 


7. 56, m 2 (2) 3 * ° 6, “(2) 


0. 288, j8 = -6. 3705, a 2 = 29. 983, and ct 2 = 10. 118. When these 

(Z) 1 ( z) z(z ) 

values are substituted into the last equation of (6) , we have u . = -68, 186. 

5(2) 

For our next approximation to r, we interpolate linearly as indicated below. 


r 


0 

-iH 

01 

i 

-624 x 10 2 

-4. 82 

-678 x 10 2 

-4. 50 

-682 x 10 2 


With r. = -4. 82, we subsequently calculate the remaining third 
( d) 

approximations just as the second approximations were calculated, except 
that this time we retain additional significant digits. We accordingly obtain 
v (3) =-23.521, p =-6.52777, m 1(3) =-7. 826, m = 3.006, 

<2,3)= 0.2775, o^ (3) 


28. 8156, o 2 . =10.3167, and finally u . , =-67, 863. 
A3) 5(3) 


Extrapolation using the values of calculated above with 


( 2 ) 


= -4. 50 


and = ~4. 82 yields r^^ = ~4. 87 as the next and final five-moment esti- 
mate of r. Corresponding to this value for r, final five-moment estimates of 
other parameters of interest become 0 * = 46.456, of = 5.352, 0* = 57.320, 
of = 3. 218 and a* = 0. 2759. 


To calculate conditional minimum chi-square estimates, we require 
values of x 2 for several values of r in the vicinity of the minimum. Since 
estimates of the basic parameters are already available for r = -4. 82 and 
r = -4. 87, we calculate expected frequencies and in turn x 2 for these values 
of r and subsequently make the same calculations for r = -4, -3, -2. 5 and -2, 
utilizing equations (29), (25), (26), (27) and (28) as previously described. 
Results of the pertinent calculations are summarized in Table IX which 
follows. 


When the values of x 2 from Table IX are plotted against corresponding 
values of r, it is readily observed that the minimum value (x 2 = 0. 72) occurs 
when r v * = -2. 65. With this value for r in ( 29) , corresponding estimates 
for the remaining parameters of interest are computed as before. Accordingly 
as final conditional minimum chi-square estimates, we find 0f * = 48. 304, 
of* = 6. 042, 0** = 57. 692, of* = 2. 951 and a** = 0. 3589. 
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TABLE IX. SUMMARY OF ESTIMATES FOR VARIOUS 
TRIAL VALUES OF r 


Parameters 



r 




Estimated 

-4. 87 

-4. 82 

-4 

-3 

-2.5 

-2 

V 

-23. 575 

-23. 521 

-22. 481 

-20.098 

-20. 001 

-19. 022 

mi 

-7. 867 

-7. 826 

-7. 146 

-6. 312 

-5. 894 

-5.475 

m 2 

2. 997 

3. 006 

3. 146 

3. 312 

3. 394 

3. 475 

p 

-6. 5535 

-6. 5278 

-6. 1344 

-5. 7287 

-5. 5644 

-5. 4317 

°-i 

28. 6397 

28. 8156 

31.7464 

35. 2933 

37. 0327 

38. 7382 


5. 352 

5. 368 

5. 635 

5. 941 

6. 086 

6. 224 

rr 2 

a 2 

10.3513 

10. 3167 

9.7791 

9. 0327 

8. 5705 

8. 0246 

°2 

3. 218 

3.212 

3. 127 

3. 006 

2. 928 

2. 833 


46. 456 

46. 497 

47. 177 

48. 011 

48. 429 

48. 848 

02 

57. 320 

57. 329 

57. 469 

57. 635 

57. 717 

57. 798 

a 

0.2759 

0. 2775 

0. 3057 

0. 3441 

0. 3654 

0. 3883 

Ms 

-67,821 

-67,863 

-68, 657 

-69,499 

-69, 835 

-70,256 








X 2 

3. 20 

3. 06 

1. 58 

0. 80 

0. 74 

0. 98 


For comparison, the expected frequencies based both on the five- 
moment estimates and the conditional minimum chi-square estimates are shown 
in Table X along with the observed frequencies. 


Agreement between observed frequencies and expected frequencies based 
on either set of estimates is satisfactory. However, in view of the large 
sampling errors inherent in the fifth sample moment, it should come as no 
surprise to find that in this example, x 2 for the five-moment estimates is 
substantially larger than that based on the conditional minimum chi-square 
estimates. For comparison, both sets of estimates calculated from the mixed 
sample are shown in Table XI along with corresponding estimates based on the 
individual components. 
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TABLE X. OBSERVED AND EXPECTED FREQUENCIES FOR 
1006 OBSERVATIONS FROM A COMPOUND NORMAL 
DISTRIBUTION 

Expected Frequencies 


Class 

Boundaries 

Observed 

Frequencies 

Based on 
Five-Moment 
Estimates 

Based on 
Cond. Min. x 2 
Estimates 

27. 5 - 31. 5 

1 ) 

’ 7 ) 

.9 j 

31. 5 - 35. 5 

5 i 

4.9 1 

5.2 ) 

35. 5 - 39. 5 

20 

21. 2 

20. 1 

39. 5 - 43. 5 

52 

53.7 

50. 8 

43. 5 - 47. 5 

86 

80. 5 

84. 6 

47. 5 - 51. 5 

98 

94. 1 

103. 3 

51. 5 - 55. 5 

200 

217. 9 

201. 6 

55. 5 - 59. 5 

363 

349. 5 

353. 8 

59. 5 - 63. 5 

164 

163. 3 

167. 8 

63. 5 - 67. 5 

15 1 

19. 5 | 

17. 5 | 

67. 5 - 71. 5 

2 ) 

.6 ) 

. 5 ) 

X 2 


3. 20 

0. 72 

d. f. 


3 

3 

P 


0. 362 

0. 868 


TABLE XI. COMPARISON OF ESTIMATES 


Parameters 

Component 

Estimates 

Moment 

Estimates 

Min. x 2 
Estimates 

e i 

47. 716 

46. 456 

48. 304 

02 

57. 607 

57. 320 

57. 692 


5.794 

5. 352 

6. 042 

°2 

3. 030 

3.218 

2. 951 

a 

0. 3320 

0. 2759 

0. 3589 



In concluding, it is deemed appropriate to emphasize that the methods 
presented in this paper are recommended only with large samples. Further- 
more, it is desirable that moments be calculated from the raw ungrouped 
data when possible. If the ungrouped data are not available, then at least 
grouping intervals should be relatively narrow in order to minimize errors 
in the higher moments from this source. When moments can only be com- 
puted from grouped data, it will usually be advisable to apply Sheppard's 
corrections. 
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ESTIMATION IN MIXTURES OF POISSON AND 
MIXTURES OF EXPONENTIAL DISTRIBUTIONS 


by A. Clifford Cohen, Jr. 

Summary 

In the analysis of experimental data, many of the distributions en- 
countered are the result of combining two or more separate component distri- 
butions. Estimation in these compound or mixed distributions is therefore of 
particular interest to aerospace scientists. Estimators are derived for the 
parameters of a compound Poisson distribution with probability density function 


-/u x -A A x 

f(x) =01 ^ — + ( 1 - a) e 


x: 


xl 


x — 0,1,2,.. . 


and for a compound exponential distribution with probability density function 
f (x) =a( i/M)e _x/M + (1 - a) (1/A)e" x/A , xiO 


where a is the proportionality factor (OSasl) and where m and X are com- 
ponent parameters. In addition to the more general case in which all parameters 
must be estimated from sample data, several special cases are considered in 
which one or more of the parameters are known in advance of sampling. 

The research reported in this paper was performed under NASA Contract 
NAS8-11175 with the Aerospace Environment Office, Aero-Astrodynamics Lab- 
oratory, Marshall Space Flight Center, Huntsville, Alabama. 


Introduction 


Many of the distributions encountered in the analysis of experimental 
data are the result of combining two or more separate component distributions. 
Accordingly, estimation in these compound or mixed distributions is of par- 
ticular interest to aerospace scientists. A previous paper [i] dealt with 
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estimation in mixtures of two Poisson distributions; these previous results 
are extended here to include several special cases wherein one or more of the 
parameters of the compound Poisson distribution are known, and in addition 
analogous estimators are derived for the parameters of the compound expo- 
nential distribution. 

The author wishes to acknowledge the assistance of Mr. Frank Clark 
for his work in establishing the IBM 7094 computer program described in the 
section entitled "Computational Procedures" and in the Appendix. 


Mixtures of Two Poisson Distributions 

The Probability Density Function. The probability density function of a 
compound distribution composed of two Poisson components with parameters 
p and X, respectively, combined in proportions a and i - a. may be written as 


-p x -X % x 

. e p . , , v e X 
f(x) = a f— + ( 1 - a) ; — 


I x = 0, 1, 2, . . . 
( OSaSl 


( 1 ) 


For convenience and without any loss of generality, we assume n > X. 

Three-Moment Estimators . The following estimating equations result 
from equating the first three factorial moments of a sample of size n to the 
corresponding theoretical moments. 


a = 


(x - X) 
(M - X) 


X0 


r = v 


[ 2 ] 


\ 


x(e 2 - r) - re = v [3] J 


( 2 ) 


where 


0 = p + X and r = jUA , 


(3) 
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and where the sample factorial moment is given by 


R 


n 


^[k] = Ti x(x- i) ... (x-k+ 1)~ , 


x=0 


( 4 ) 


in which R is the largest observed (sample) value of the random variable x, 
n is the sample frequency of x, and 


R 



For simplicity of notation, x has been written in place of for ti 16 first 
sample factorial moment. 

On solving the last two equations of (2) simultaneously for r and 9, 
it follows that 

\ 


(5) 


where the asterisk (*) distinguishes estimators from the parameters being 
estimated. The required estimators of p and X follow as 


0 * - " [ - 3] ~ X > ] 
V [2 ] " X 


r* = x6* - v 


[ 2 ] 


p* = ~ [ 0* + \iW - 4r 5 " ] 

Li 


\ 


a* = | [ 0 * - N le* 2 - 4r* ] j 


( 6 ) 
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These estimators are the two roots r t and r 2 of the quadratic equation 


Y 2 - 0*Y + T* = 0 , ( 7) 


where p* = and A* = r 2 , (rj > r 2 ). The proportionality parameter a is 
estimated from the first equation of ( 2) as a* = (x - A' 1 ' )/ ( M' 1 ' - A* ) . 

The estimators given in equation (6) were originally derived by 
Rider [13] , but he employed ordinary rather than factorial moments with the 
result that his derivations were somewhat complicated and his expressions for 
0* and were more involved than those given here. 

Estimators Based o n the First Two Sample Moments and the Sample 
Zero-Frequency . It is well known that the higher sample moments are subject 
to appreciable sampling error, and in an effort to reduce errors from this 
source, the estimating equation based on the first two sample moments and 
the sample zero-frequency, was derived [14] as 


x - A nty/n - e 

G( A) - A -G( A) -A 

0 “ G 


( 8 ) 


in which 


G( A) = 



(9) 


where n 0 is the sample zero -frequency. Equation ( 8) can be solved for 
A** using standard iterative procedures and, with A* v thus determined, esti- 
mators of p and a follow as 


,* * 


[ 2 ] 


— -v 

- xA 


x - A"* 


a 


*>:< 


x - A v * 

jLt** - A** 


( 10 ) 
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The double asterisk (**) distinguishes these estimators from the three- 
moment estimators and in turn from the parameters being estimated. 
Unfortunately, no simple procedure for solving equation (8) has been devised. 
However, a computer program based on iterative procedures described by 
Whittaker and Robinson [ 15 (Chap. VI) ] has been developed (see Appendix) 
to solve equation (8) , using as a first approximation the three-moment esti- 
mate of A given by equation ( 6) . 

Estimation With Some P arameters Specified . 

a. a Known 


In this case, we need only estimate m and A; for this purpose, 
the first two equations of ( 2) may be written as 


x - A \ 



x(p + A) - pA = v j 


( 11 ) 


where Q and r have been replaced by their defining relations as given in 
equation (3). 

With a known, we obtain the following quadratic equation in A from the 
two equations of (11): 


x 2 - av 


>2 _ OT7 


2xA + 


[ 2 ] 


1 - a 


= 0 . 


( 12 ) 


On solving equation (12) 



(13) 
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and from the first equation of ( ii) 


/J* = [x - X* ( 1 - a.) ]/a . ( 14) 

b. a and p Known 

In this case, X may be estimated from the first equation of (11) as 


_ x - ajx 


1 - a 


(15) 


c. a and X Known 


In this case, it follows from equation ( 11) that 


* = x - (1 - oQX 
^ a 


(16) 


d. m Known 

In this case, we may employ equations (11) to estimate a and X. 
Accordingly, from the second equation of (11) 


x*-4!l 


- xM 


X - jU 


(17) 


and from the first equation of (11) 


* x - X* 
a - 


M - X* 


e. X Known 

In this case, the second equation of (11) gives 


(18) 
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- xA 


( 19 ) 


M* = 


y [2] 

x - A 


> 


and from the first equation of (11) 


Mixtures of Two Exponential Distributions 


( 20 ) 


Th e Probability Density Function. In many respects the exponential 
distribution may be thought of as a continuous analog to the discrete Poisson 
distribution. In any event, estimating equations in mixtures of two exponential 
distributions quite closely parallel the estimating equations considered in the 
preceding section for mixtures of two Poisson distributions. Consider a com- 
pound exponential distribution with probability density function 


f(x) = a(l/ju)e X/M + (1 - oi) ( 1/A) e x/X 


x ^ 0 

< n > X > 0 
0 S a S 1 


( 21 ) 


The nonessential restriction that /i > A is imposed as a matter of convenience 
and without any loss of generality. 

The kth noncentral moment of x is 

OO 

m' = f x k f (x) dx = k! [afx*’ + ( 1 - a) A^] . (22) 

k 0 

Accordingly, the first three noncentral moments are 

m| = ap + ( 1 - a) A 
m] = 2[ct^ 2 + (1 - a) A 2 ] 
mi = 6[an 3 + ( 1 - a) A 3 ] 
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Three-Moment Estim ators. When the first three noncentral sample 
moments, designated p[, and v%, respectively, with v\ = x, are equated to 
the theoretical moments of ( 23) , we obtain the estimating equations 


x - A = a (fj, - X) 


f- - A 2 = a( » 2 - A 2 ) 


g - a 3 = o?(p 3 - a 3 ) j 


(24) 


These equations differ from the corresponding equations for mixed Poisson 
distributions only in that v{/2 and v' 3 /6 have replaced the factorial moments 


[ 2 ] 


and 


of the mixed Poisson distribution. 


On eliminating a between the first and second and between the first and 
third equations of ( 24) , we simplify to obtain 


xO- r = 



x(0 2 - r) - re 



(25) 


which are completely analogous to the last two equations of ( 2) in the case of 
mixed Poisson distributions. Here, as in the Poisson case, 0 and T are defined 
by equation (3). Accordingly, on solving the two equations of (25) simul- 
taneously, we have as estimators of 0 and r 



/ 


( 26 ) 


which are analogous to equation (5) for the mixed Poisson distribution. 
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Finally, with 0* and r* determined from (26), n* and X* follow 
from equation ( 6) as in the Poisson case , and a* follows from the first 
equation of ( 24) as 


a 


* = 


-x * 

x - A 

At* - X* * 


( 27 ) 


Estimation With Some Parameters Specified. 
a. a Known 


We need only replace v ^] " 2/2 an< * quadratic equation of 

(12) becomes, for the present case, 


X 2 - 2xX + 



0 . 


Accordingly, 



+ x - X* ( 1 - a ) 



(28) 


(29) 


b. a and m Known 

In this case, the estimator for X follows from the first equation of 

(24) as 


x * = x z_ay+ 

1 - a ’ 

which is identical with the corresponding estimator, equation (15) , in the 
Poisson case. 


(30) 
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c. Q? and A Known 


In this case, it follows from the first equation of (24) that 


x - ( 1 — a) A 

a 


( 31 ) 


d. A Known 


In this case, we need only replace ln e( l ua ^ on (19) with 

1 ^ 2/2 and, accordingly, 


1* = 


}± 

_2 

x - A 


xA 


(33) 


a* = 


x - A* 
M - A* 


Computational Procedures 


The solution of the transcendental estimating equation ( 8) from the 
section entitled "Mixtures of Two Poisson Distributions" provides an interesting 
illustration of iterative numerical computational techniques described by 
Whittaker and Robinson [ 15] . To facilitate solution of equation (8) , the de- 
nominator of the left side is interchanged with the numerator of the right side, 
and the resulting equation becomes 


x - 


n 0 /n - 



e 


G(A) - A 

-G( A) -A 

- e 


(34) 


where G(A) remains as given by equation (9) . 
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Equation (34) might be condensed to the form L(X) = R(X) where 


L(X) = — X ~ and 

n 0 A - e 


R(X) = 


e 


G(X) - X 
-G(X) „ 


-x • 


( 35 ) 


The two functions L(X) and R(X) are essentially as plotted in Figure 1 

below. 


L(A) 



We begin with an initial approximation X 0 and iterate toward the value 
X** as described by Whittaker and Robinson [15 (pp. 81-83) ]. The three- 
moment estimate of X given by equation ( 6) of the section entitled "Mixtures 
of Two Poisson Distributions" provides a satisfactory value for X 0 . This 
initial approximation is substituted into the second equation of (35) to obtain 
R 0 , which is merely an abbreviated notation for R( X 0 ) . We then solve the 
equation 


L(Xi) = R 0 


(36) 
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to obtain A 1( the next approximation. This cycle is repeated as many times as 
necessary to attain the desired degree of accuracy. Equation (36) is itself a 
transcendental equation, though somewhat simpler in form than the original 
equation ( 34) . It is amenable to solution by the Newton-Raphson method 
[15 (pp. 84-86) ] . For the ith cycle of iteration, the equation corresponding 
to (36) becomes 


L( A.) 


x - A. 

l 

n(/n - e ^ 


R 


i-1 


(37) 


which may be written as 


f(A.) = 0 , 

where 

f( A.) = A. - R. , e" Ai - C. (38) 

i l l-l l-l 

and 

C._ 1 = (x - R._ 1 n,/n) . 

Equation (37) may be readily solved using the Newton-Raphson method, 

where A. . , the ( r + 1) st iterant to A. , is given by 
i:r+i i 

f( A. ) 

1 = A — . 

i:r+l i:r f' (A. ) 

i:r 

The first derivative of f(Aj) follows from equation (38) as 


f'Uj) = 1 +R._ 1 e 


-Ai 
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Accordingly, 


\:r+l \:r 


A. - R. . e 
i:r l-l 


_ ^i:r _ 


C. 


i -1 


1 + R. , e Xi 
l-l 


(39) 


As an initial approximation A. to A^, it will usually be satisfactory to 

let A. ^ = A The Newton-Raphson iterative technique is continued through 

as many cycles as necessary to attain the desired accuracy in A.. More 

specifically, this procedure is terminated at the end of the rth cycle, the 
first cycle for which 


where specifies the maximum permissible absolute value deviation. With 
A^ thus determined, we calculate R^, set up the new equation 


L( A. ,) = R. 
l+l l 


and continue the primary routine through k cycles. The kth cycle is the first 
for which 


|L k - R R | < 63 , (40) 

where 6 2 specifies the maximum allowable absolute value deviation. The re- 
quired estimate of A is then 



Illustrative Examples 


Mixed Poisson Distribution . To illustrate the application of his three- 
moment estimators, Rider [13] chose an example constructed by mixing equal 
proportions of two Poisson distributions with M = 1. 5 and X = 0. 5, respectively. 
These data are as follows: 


X 

0 

1 

2 

3 

4 

5 

6 

7 

n 

X 

830 

638 

327 

137 

49 

15 

3 

1 


In summary, n = 2000, n 0 = 830, x = 0. 9995, v ^ j = 1. 248 and = 1. 734. 

Direct substitution of these values into equations ( 5) and ( 6) yields the three- 
moment estimates 


M* = 1.4766563, 

X* = 0. 47765894, 
a* = 0. 52236479. 

The above results differ slightly from those given by Rider caused apparently 
by small round-off errors in his calculations. 

Estimates based on the first two sample moments and the sample 
zero -frequency, calculated by a computer program of the routine described 
in tiie preceding section, are 


/Li** = 1. 4936, 

X** = 0. 4956, and 
a** = 0. 5049. 
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These estimates are in much closer agreement with the actual population 
parameters p = 1. 5, X = 0. 5, and a = 0. 5 than the three-moment estimates. 
Investigations are continuing with regard to the relative efficiency of the 
three-moment and the two-moment plus zero-frequency estimates; but at 
least in the present instance, where a large proportion of the population is in 
the zero class, the two-moment plus zero-frequency estimates seem to be 
preferred. 

Mixed Exponential Distribution . To illustrate the application of esti- 
mators derived in this case, a sample of 2000 observations was selected from 
a mixed population constructed by combining two exponential distributions with 
H = 2, X = 1, and a = 0. 4. Data for the sample selected are summarized as 
follows: n = 2000, x = 1. 42, = 4. 38, and = 21. 6. 

Direct substitution of these data into equations (26) , (6) , and (27) 
yields as three-moment estimates: 


M* = 1. 85, 
X* = 1. 02, 
a* = 0. 48. 
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Appendix 

FIND - A Computer Program* 

FIND is a Fortran IV computer program which calculates estimates for 
the parameters a, /it, and A of a compound (mixed) Poisson distribution. 
These estimates are calculated from (a) the first three sample moments and 
(b) the first two sample moments and the sample zero-frequency. 

In finding X for the second case, the following equation is solved: 


x - X nn/n - e 

G(A) - X -G(A) -X ’ 
e - e 


where 


G( A) = 


'[ 21 ^ 


x - X 


and v ^] ’ x ’ an< * n o/ n are known constants. FIND makes use of the Newton- 
Raphson and geometrical iteration methods [15] in solving the equation. 

FIND requires, for each data sample, input values for x, n</n, > 

and , punched on a single card. Iteration continues through k < 500 

cycles until the absolute error of equation (40) is less than 0. 00001, i. e. , 
until 


I L k - R k I < 0. 00001 . 

If this criterion is not met when k = 500, the message "completed 500 itera- 
tions with no success" is given and the program stops. Should greater 
accuracy be required in the estimate of X, appropriate change should be 
made in the source program card "TOL = 0. 00 . . . . " 


* The computer program presented was by Mr. Frank C. Clark. 
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FIND prints out the following: 


1. Values of the index, i. 

2. Values of A in the Newton-Raphson iteration. 

3. Values of 

ERROR = TEST 1 - TOL, 
where 

TEST 1 = |L. - R. . I. 

i:r i-1 

4. a, n, and A based on the first three sample moments. (This 
value of A is used as the first approximation in the Newton-Raphson process. ) 

5. a , p, and A based on the first two sample moments and the 
sample zero-frequency. 
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FIND (FORTRAN IV) 


ESTIMATION IN MIXTURES OF TWO POISSON DISTRIBUTIONS 

DIMENSION LAM(AOOO) 

REAL MU. NUE2 > NUE3. NON. L AM. L. LAMBDA. MU1 

1 READ(5*2)XBAR.NON*NUE2*NUE3 

2 FORMAT ( 4F 10 . 5 ) 

THET = <NUE3-XBAR*NUE2) / ( NUF2- ( XBAR**2 ) ) 

CLAM = XBAR*THET-NUE2 

MU = (THET+SQRT(THET**2-4.0*CLAM) )/?.0 
1=1 

LAM ( I)= (THET-SORT (THET**2-4.0*CLAM) )/ 2.0 
N = 0 

ALPHA 1 = ( XBAR-L AM ( I ) ) / ( MU-LAM ( I ) ) 

K = 0 

G = <NUE2-XBAR*LAM< I ) ) / (XBAR-LAM(I) I 
N = N + 1 

R = { G— L AM ( I ) ) / ( EXP ( — G ) — FXP ( —LAM ( I ) ) ) 

9 C = ( XBAR-NON*R ) 

10 K=K+1 

LAM ( 1 + 1 ) = LAM ( I ) - ( ( LAM ( I ) -R*EXP ( -LAM ( I ) ) -C ) / ( 1 . 0+R*FXP ( -LAM ( I ) ) ) ) 
L = ( XBAR-L AM { I +1 ) ) / ( NON-EXP (-LAM( 1+1 ) ) ) 

TOL = .00001 
TEST 1 = ABS(L-R) 

C 

60 FORMAT ( 1H . I 5 . 5 X » El 5 . 8 . E 1 5 . 8 ) 

ERROR = TEST 1 TOL 
WRITE (6. 60) I . LAM ( I ) .ERROR 
IF (TEST1-TOD20.15.15 

23 1=1+1 

GO TO 10’ 

20 G = ( NUE2-XBAR*LAM ( 1+1) )/(XBAR-LAM( 1+1) ) 

R = (G-LAM( 1+1) )/(EXP(-G)-EXP(-LAM( 1+1 ) ) ) 

TEST? = A.BS ( L-R ) 

IF (TEST2-TOL130. 25*25 

24 1=1+1 
K = 0 

GO TO 9 

C 

15 IF( 500-022*22*23 

25 IF ( 500-NJ22 *22.24 
22 WRITEI6.28) 

28 FORMAT < 42H1C0MPLETED 500 ITERATIONS WITH NO SUCCESS) 

GO TO 100 

30 MU1 = ( NUE2-XBAR*LAM( 1+1 ) ) /(XBAR-LAM( 1+1 ) ) 

LAMBDA = LAM ( 1+1 ) 

ALPHA2 = (XBAR-LAMII+l) )/ CMU1 -LAMII+1)) 

C 

WR I TE ( 6 . 50 ) 

50 FORMAT (39H1ESTIMATES BASED CN FIRST THREE MOMENTS) 
WRITE(6.51)MU.LAM( 1 ) , ALPHA 1 

51 FORMAT ( 10H0 MU =E15.8»10H LAMBDA = E15.8.9H ALPHA = E15.8) 
WRITE(6.52) 

52 FORMAT (74H0ESTIMATES BASED ON FIRST TWO SAMPLE MOMENTS AND THE ZER 
10 SAMPLE FREOUENCY) 

WRITEI6.53) MU1 .LAMBDA . ALPHA2 

53 FORMAT ( 1 1H0 MU = E15.8.10H LAMBDA = EJ5.8.10H ALPHA = E15.8) 

GO TO 1 

100 STOP 
END 
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ESTIMATION OF PARAMETERS IN COMPOUND 
WEIBULL DISTRIBUTIONS 

by Lee W. Falls 

Summary 

The two-parameter Weibull distribution has been recognized as a useful 
model for survival populations associated with reliability studies and life testing 
experiments. In the analysis of atmospheric data, the distributions encountered 
are often a result of combining two or more component distributions. These 
compound distributions are consequently of interest to aerospace scientists. 

This paper presents a method for estimation of the parameters of a compound 
Weibull distribution with density function. 


f(x) = afj(x) + (i - a) f 2 (x) , 0 < a < 1 


where 


f t (x) = Ti 1 x^ 1 * exp [-x Tl /0i] , x ^ 0, > 0, y l > 0 


and 


f~l 


The parameters required are a, the proportionality factor, y it y 2 , 9 j and 0 2 . 
The most general case of estimation will be considered in addition to a number 
of special cases that may be of practical value. 

Introduction 

The Weibull distribution, derived in 1939 by W. Weibull, has been 
recognized as an appropriate model in reliability studies and life testing. 
Numerous methods for obtaining efficient estimates of the two parameters 
of this distribution have been outlined in recent years [16-18, 20] . 

In actual physical applications, however, a mixture of two Weibull 
distributions often seems to be a more desirable model. Distributions re- 
sulting from mixing two or more component distributions are designated as 
"mixed" or "compound. " This situation is quite common in the analysis of 
atmospheric data and consequently is of interest to aerospace scientists. 
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Compound normal, Poisson and exponential distributions have been studied 
by A. C. Cohen, Jr. [2, 3,19]. A method for estimating parameters of mixed 
distributions using sample moments has been outlined by Paul R. Rider [13] 
who considered compound Poisson, binomial, and a special case of the com- 
pound Weibull distribution. A graphical procedure for estimation of mixed 
Weibull parameters in life-testing of electron tubes is given by John H. K. 
Kao [20] . Although graphical methods have value for locating outliers, 
deriving initial estimates, and for determining whether the distribution is as 
hypothesized, for estimation purposes the analytic approach is probably 
superior. 

This paper represents an attempt at estimating, by the method of 
sample moments, the five parameters of the compound Weibull distribution 
with density function 


f(x) = cyf x) + (1 - a) f 2 (x) , 0 < a. < 1 (1) 

where 

fj(x) = 7*01 I x^ 1 * exp [-x^V^i ] » x ^ 0, 0* > 0, 

f2< x ) = 72^2 lx?2 1 ex P [-x Y2 /0 2 ] , x ^ 0, 0 2 > 0, 

The parameters involved are two scale parameters 0j and 0 2 , two shape 
parameters yj and y 2 , and the proportionality parameter a which expresses 
the probability that a given observation x. comes from the population fj. 

The compound cumulative distribution function is defined 

F(x) = aF^x) + (1 - a) F 2 (x) = i - a exp [-x^V^il 
- (1 - a) exp [-x^ 2 /0 2 ] • 

Figures 2 and 3 illustrate a generalized mixed Weibull probability 
density function and its corresponding distribution function. 



Tl > 0 

72 > 0 


( 2 ) 
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The most general case of estimation will be considered in which all 
five parameters must be estimated from the data. Also, a number of special 
cases will be investigated in which certain parameters are known in advance 
of sampling or are restricted in some manner. Included will be the special 
case where y* = i, i. e. , 

fj(x) = 0J 1 exp [-x/ej , (4) 


which is the well-known exponential distribution. Estimating procedures are 
greatly simplified in these special cases as there are fewer sample moments 
involved in the estimating equations. Also, sampling errors are reduced 
because of the elimination of the need for higher order moments. 


Estimation in the General Case 

The rth theoretical moment about the origin of f(x) is given by 


OO OO 

ju '= a J x r fj(x) dx + (1 - a) J x° f 2 (x) dx , (5) 

r 0 0 


where fi(x) and f 2 (x) are defined as in equations (2). The first five theo- 
retical moments about the origin of ( 1) follow as 


>*1 = 0 e l /T>r (^" + l) 

A - « of /n r (^ + 


2 1/72 r ( 

— + 1 ) 


A \ 

-72 ' 


l? /y * T 

(- +1 ) 





i 2 3/ ^ r 

(- +1 ) 



\72 / 

/ 


(Equation (6) concluded on following page. ) 


( 6 ) 
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Mi = 

a r 1 

(n + 1 ) 

+ (1 - a) et /72 r | 

(M 

M 5 = 

0, e\ /y ' r 


t + (i - «) 0 2 5/Y2 r j 



I 

> ( 6 ) 
(Con- 
cluded) 

/ 


where r is the gamma function, i.e. , 


OO 


r<k) = / 

0 


k-i 

y 


e" y dy . 


Employing the technique of equating population moments to correspond- 
ing sample moments, the set of equations (6) becomes 


mj = a 0^ Yl T i) + (1 - a) O^ 72 T ^~+ i) 

mi = a e 2 / 7i r (^- + i) + (i - a) ef /72 r (^-+ 1 ) 

m' 3 = a e 3 / 71 r(^-+ l) + (1 -a) e 3 / 72 r (^+ l) 

m i = a S? /Yl r + i) + ( 1 - a) ef /72 T + l) 

m 5 = Q! sf /Yl r (~ + l) + <i-«> e 2 5/Y2 r + i) 


> (7) 


where mj (i = 1, 2 5) is the ilh noncentral moment of the sample. 
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The set of equations (7) is a system of five equations which must be 
solved simultaneously for estimates of the five parameters a, 0 it 0 2 , yj and 
y 2 . For convenience in handling these equations, we will make the following 
transformations where necessary in this paper. (This notation will be used 
unless stated otherwise. ) 

Let 


- e\ /yi 

Pi = r 


^i = r, 

(M 

> 


II 

to ^ 

Pi = r 

(M 

= r 

(,T +1 ) 


► (8) 


P 3 = r 


^3 = r 1 

GM 




Thus, the first three equations of (7) become 


mi = av| 3 t + (i - a) nipt 
m 2 = av 2 /3 2 + (1 - a) n 2 ifc 
m 3 = av 3 / 2 3 + ( 1 - a) n 3 ip 3 . 


(9) 

( 10 ) 
(ID 


Solving (9) for v, substituting this expression into (10) , and then solving 
for n yields 


mj( i-a ) i/)j/3 2 ± n /- ( m i) 2 o! ( i-a ) /3 2 /? 2 ^ 2 + m^o; ( i-a) 2 (3 2 p\ip\ + m 2 Q' 2 ( i-a)(3\jp 2 

n = — — — 

( l-a) 2 i/ij/3 2 + a ( 1-a) fS\ip 2 


( 12 ) 
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Substituting the expression for n from ( 12) back into ( 9) and solving for 
v gives 


v = ■ 


mj - (i-oi)^i mj(i- 

G^i 


a) ip-fez ± *s/- ( m J) 2 a ( 1-a ) + mga( i-a) 2 ^\>Pi + 


(i-a) 2 ^ + a(i-a)/3^ 2 




(13) 


Upon substituting the expression for n from ( 12 ) and the expression for v 
from ( 13) into equation ( 11 ) , we have one equation in the three unknowns 
a, and y 2 . At this point, it is obvious that explicit expressions for the un- 
known parameters a, yj and y 2 cannot be obtained. Therefore, it is suggested 
that the following procedure be used to obtain a graphical estimate of a , the 
proportionality parameter. This method is essentially that of Kao [ 20 ] and is 
based upon the fact that a simple Weibull cumulative distribution becomes a 
straight line in In versus ln-ln coordinates. This method of estimating a 
will produce a relatively small error in the estimating procedure since a is 
limited to the range 0 < a < 1. 


a. Plot the sample cumulative distribution function for the mixed data 
on special In versus ln-ln paper and visually fit a curve (called Weibull plot) 
among these points. * 


b. Starting at each end of the Weibull plot, draw two tangent lines and 

denote them by QfFi and (1 - a)F 2 which are estimates of o;Fi(x) and 
( i - a) F 2 (x) , respectively. 


c. At the intersection of ( 1 - a) F 2 with the upper borderline drop 


a vertical line whose intersection 
gives our estimate of a. 


with aFj as read from the percent scale 


See Figure 5 under the section entitled "An Illustrative Example" for an 
illustration of this method. 


Once an estimate of a has been determined graphically, solve equa- 
tion (11) for y t and y 2 by the following iterative procedure. This procedure 
is taken from Cohen [3] and is a modified Newton-Raphson method. 


* Special Weibull graph paper is available from Cornell University, Ithaca, 
New York. 
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Assume a value for and solve equation (11) for a first approxi- 
mation to y 2 . These first approximations can be substituted into ( 12) and ( 13) 
to obtain first approximations to 0j and 0 2 . The first set of approximations 
is then introduced into the fourth equation of ( 6) to approximate the fourth 
noncentral theoretical moment, 


Let denote the ith approximation to y t and let p^ denote the 

ith approximation (corresponding to y^.^) to pj. It should be relatively 
easy to find approximations y^^ and y^ + ^ such that the sample moment 
mj is in the interval [p^.^ , p^ , + ^]. Once the interval between y^^^ and 
y . has been narrowed sufficiently, the required estimate y^ can be ob- 

-L V l - * J - } 

tained by a simple linear interpolation as indicated below. 


7i 

Ml 

T l(i) 

M 4(i) 

n 

m 4 

y i(i+D 

M 4(i+1) 

J 


The required estimate of y 2 can subsequently be obtained from equation (11). 
Once Yi and y 2 have been determined by equation (11) , estimates for and 
0 2 are obtained from equations (13) and (12), respectively. 

Unfortunately, the quadratic solutions in equations (12) and (13) 
result in more than one set of estimates. The problem of non-unique sets 
of estimates was considered by Karl Pearson [4] and A. C. Cohen, Jr. [3] 
in connection with mixtures of two normal distributions. Pearson suggested 
choosing the set of estimates which gives closest agreement between the sixth 
sample moment and the sixth theoretical moment after equating the first five 
sample moments to the corresponding theoretical moments. This procedure 
is followed for all acceptable sets of estimates. 

Cohen [3] suggests, as an alternate procedure for resolving the 
problem of multiple sets of estimates, that we might choose the set of esti- 
mates which produces the smallest Chi-Square Index of Dispersion when 
observed frequencies are compared with expected frequencies. 
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In the general case of estimation considered here, we are concerned 
with a mixture of two Weibull distributions where the proportionality factor a 
is estimated graphically from the cumulative frequencies and the four re- 
maining parameters are estimated by equating the first four sample moments 
to corresponding theoretical moments. When confronted with more than one 
set of acceptable estimates, we adopt Pearson's suggested procedure and 
choose the set which produces the closest agreement between the fifth non- 
central moment of the sample m§ and the theoretical moment M 5 given by 
the final equation of ( 6 ) . 

The calculations described above may be carried out iteratively with 
relative ease using the computer program included in this paper as an 
Appendix. First approximations to initiate the iterative process may be ob- 
tained using a graphical method such as that of Kao [20.]. 

Estimation in Special Cases 

A number of special cases that may be of practical value in which 
certain parameters are known or are restricted in some manner are con- 
sidered. 

fli Known . With 9i known, we must estimate the parameters ot, 0 2 » 

1 /vi 

and y 2 only. If we let v = 6\ , equations (9), (10) and (11) become 

m] = a. 9Y' yi p i + (1 - a)wp t , (14) 

m 2 = a e^ Yi p 2 + U - ot) n 2 ^ , (15) 

m 3 = a O^ 71 p 3 + (1 - a) n z 4> z . (16) 


Solving (14) for a gives 


_ mj - ni/>i 

dY^Pi- nft 
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Inserting this expression for a into equation ( 15) and solving for n, we obtain 
after considerable algebraic manipulation 


- m-^, ± n/IhiJ) 2 ^ 2 ! - 2m 2 v 2 ifi/3 2 + vVt^l - 4mlv 3 (l l /3 2 >h * + 4m|v 2 (3 t 2 ifc - 4mjmJv/3 1 if 2 

2v/3 t fe - 2mi(fe 


(18) 


Substituting this expression for n back into ( 17) gives 

|~v 2 i 


m\ - ipi 


v 2 ^ - m zh ± ^ (ni 2 ) 2 ^i - 2m2vV 1 /3 2 + v 4 $/S§ - 4m\v^^ 2 h + 4<mi) + 4m^v 2 /3 2 1 ^ 2 - 4m}m 


«1 1/Yl ^i - - 


2v0ifc - 2m ^ 


v 2 ^ 2 " m 2$i i's/tmj) 2 # - 2mlv 2 $p2 + v 4 ^l - iniv 3 /3^ 2 ^ + + ^mjv 2 /^ - Wpnjv/S^ (19) 

ZvPtfk - 2m t i ip 2 


Now, upon substituting the expression for n from (18) and the expression for 
a from (19) into equation (16) , we have an equation in the two unknowns yj 
and y 2 which may be solved by the iterative procedure described in the general 
case. Once y* and y 2 have been determined, we obtain our estimates for 0 2 
and a from equations (18) and (19) , respectively. As in the general case, 
the positive and negative roots resulting from the quadratic solution in (18) 
unfortunately results in more than one set of estimates. As before, it is 
suggested that the set of estimates which gives the closest agreement between 
the fifth noncentral moment of the sample and the corresponding "fitted" com- 
pound curve be used. 

6 9 Known . With 0 2 known we need only estimate a, 0 lf y 1 and y 2 . 

If we let n = 0 2 ^ 2 , equations (9) , (10) and ( 11) become 

m[ = avfi t + ( 1 - a) 0 2 /y2 , (20) 

m^ = av 2 |3 2 + ( 1 - a) 0% // ' Y2 , (21) 

m 3 = Q?v 3 /3 3 + (1 - a) 0 2 / ' y ‘ 2 ip 3 . ( 22 ) 
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Solving (20) for a gives 


ml - 0^21 h 
v0i - oy yz ipi 


(23) 


Inserting this expression for a into equation (21) and solving for v, we obtain 
after considerable algebraic manipulation 


m^i - n 2 ( 3 $2 ± ^n‘/3 - 2n 2 mJ/3 2 ,^ 2 + (mj) 2 /^ - + + 4(m' 1 ) 2 aV ! ifc - 4m 

V = — 

2 mJ/ 3 2 “ ZntyiPz 


(24) 


Substituting this expression for v back into (23) gives 


mj - 


r mjpi - n 2 #,^ ± N/n 4 ^!/! - 2n J m2/3 2 & + (m 2 ) 2 /?| - 4m;in 2 n tpfa + + 4(m'i) 2 n ! p 2 i/^ - 4m' 1 n 3 !f>i&0 2 


*1/ 


2mJ/3 2 “ 


- - 9. 


(25) 


'2 Vl 


Now, upon substituting the expression for v from (24) and the expression for 
a from (25) into equation (22) , we obtain an equation in the two unknowns y j 
and y 2 which may be solved by the iterative process described previously. 

Once yj and y 2 are determined, we obtain our estimates for and a from 
equations (24) and (25). As in the case of known, the positive and nega- 
tive roots resulting from the quadratic solution in (24) gives more than one set 
of estimates. Again, it is suggested that we choose the set of estimates which 
gives the closest agreement between the fifth noncentral moment of the sample 
and the corresponding n fitted M compound curve. 

yi Known . If yj is known we must estimate a, 0j, 0 2 and y 2 only. 
Solving equation (9) for v gives 


v 


m} - (1 - a) n^i 
<*Pi 


(26) 
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Inserting this expression for v into equation ( 10) and solving for n, we 
obtain 


m i( i-a ) ij> j/3 2 ± \l - ( m J) 2 a ( i-a ) P 2 p\rp 2 + m 2 a ( i-a ) 2 p 2 p\ip\ + m 2 a 2 ( 1-a ) pfip 2 

n = ... — “ .. — - — — . 

(i-a) 2 $3 2 + &(l-oi)p\ip 2 

(27) 

Substituting this expression for n back into (26) gives 


mj - ( t-a)ip 1 

v = 

a/3j 


m}( l-a)4>jP 2 + m^U-aO/s}^ 

(l-a) 2 ify3 2 + a( l-a)/3^ 2 




(28) 


Upon substituting the expression for n from (27) and the expression for v 
from (28) into equation ( 11) , we obtain an equation in the two unknowns a and 
y 2 which may be solved by the iterative procedure described in the general 
case. With a and y 2 determined we may solve equations (27) and (28) for 
0 2 and 0 lf respectively. As before, the positive and negative roots which 
result from the quadratic solution in (27) give more than one set of estimates. 
Again, we choose the set of estimates which gives the closest agreement to 
the fifth noncentral moment of the sample. 

An alternate method for estimation in this case would be to estimate 
a. graphically as in the general case and then solve equation (11) for y 2 after 
the substitution of the expression for n from (27) and the expression for v 
from (28) into equation (11). As before, 0 2 and d 1 would then be obtained 
from equations (27) and (28). 

y? Known . With y 2 known, we must estimate a, 0 1# 0 2 a nd y\ only. 

As in the case, y\ known, solving equation (9) for v gives 


m} - (1-a) wpx 

a Pi 


(29) 
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Inserting this expression for v into equation ( 10) and solving for n, we 
obtain 


mj( 1—ot) 1P1P2 ± \/ - ( m }) 2 a ( 1 -a ) (3 2 p\ip 2 + m 2 o! ( 1 -a ) 2 P 2 p 2 ip\ + m 2 a 2 ( 1- a ) (3\ip 2 

n = 

(l-a) 2 $/3 2 + a(i-a)p\ip 2 


(30) 


When we substitute this expression for n back into (29) , we get 


mj - ( l-a)ip i m}( 1 

oil3 1 


-a)tp^ 2 ± ^ -(m\) 2 a(i-a)/3 2 ^\>h + mj«( l-a) 2 p^\ip\ + 


( l-a) 2 i/ijP 2 + a(i -a)p\ip2 




(31) 


Upon inserting the expression for n from ( 30) and the expression for v from 
(31) into equation ( 11) , we obtain an equation in the two unknowns a and y* 
which may be solved by the iterative procedure described previously. With 
a and yj determined in this manner, we now solve equations (30) and (31) 
for 0 2 and respectively. As before, we choose the set of estimates 
which gives the closest agreement between the fifth noncentral moment of the 
sample and the corresponding "fitted" compound curve. 

As in the case of yj known, an alternate method for estimation would 
be to estimate a, the proportionality parameter graphically, and then solve 
equation ( 11) directly for y* after the substitution of the expression for n 
from (30) and the expression for v from (31) into equation (11). 

Ys = 1 . This is a special case of yj known. Thus, the case is reduced 
to mixing an exponential distribution with a Weibull distribution where fj(x) 
in equations (2) is an exponential distribution and falx) is a Weibull distri- 
bution. We need only estimate a , y 2 , 0j and 0 2 . 

With yj= 1, equations (9), (10) and (11) become 


(32) 

(33) 


mj = av + ( 1 - a) mpi 
m 2 = 2 o!v 2 + ( 1 - a) n 2 ip 2 
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m3 = 6av 3 + ( i - a) a 3 fa , 


(34) 


where the only change in notation from previous cases is v = 
Solving equation ( 32) for v gives 

m } - ( 1 - ol) nfa 


(35) 


Substituting this expression for v into equation ( 33) and solving for n gives 


2mJ(i -ot)ipi ±\l2m20i(i-a) 2 4)\ - 2(m[) 2 QL(l-cL)fa + m^a 2 * 1 -ot)ip 2 

n = — . (36) 

2 ( 1 -a) 2 ip\ + ol ( 1 -a) fa 


Inserting this expression for n back into (35) gives 


v = - i Pi 

a 


2m J( l-o:)ip 1 ± \j 2 m 2 ,a(i-a) 2 ip\ - 2(mJ) 2 Q!( i-o>) fa + ni2G! 2 ( 1-a) fa 


2a( i-a) ip\ + a 2 fa 


(37) 


Now, upon substituting the expression for n from (36) and the expression for 
v from (37) into equation (34) , we obtain an equation in the two unknowns ol 
and y 2 which we may solve using the iterative procedure described in previous 
cases. With minor simplifications, equation (34) becomes 


mj = 6a 


mi 

~ $1 

a 


(\ 2mJ( i-a) fa ± n/ 2 m 2 a( i-a) 2 ^ - 2(mJ) 2 a( i-a) fa + m 2 a 2 ( i-a) fa 
^ 2a( i-a) + a 2 !^ 



( 38 ) 


(i-a) fa 


2mJ(i-a) ipi ± \l2mlai i-a) 2 ip\ - 2(mJ) 2 a ( i-a) fa + rnia 2 (i-a)fa 


2( i-a) 2 ^j + a( i-a) fa 
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Once a. and y 2 have been determined from equation (38) we may obtain 
our estimates for 0 2 a nd from equations (36) and (37) , respectively. 

As in the other cases, we choose the set of estimates which gives the 
closest agreement between the fifth noncentral moment of the sample and the 
corresponding "fitted" compound curve. 

An alternate method of estimation would be to estimate a graphically 
as in the general case and then solve equation (38) directly for y 2 . As above, 
estimates for 0 2 and 6 j would then be obtained from equations (36) and (37). 

Yi 72 ~ Unknown . Changing our notation slightly, we will let 
7i = y 2 = y. Thus, we must estimate a, y, 0j and 0 2 . Now, let 


v = f>K T 

II 

•n 

(H 


'r-f CSJ 

II 

S3 

/3 2 = r | 

(H- 



With this notation, equations (9), (10) and (11) become 


mi = ayfii + (i - a) n/^ 

(39) 

m 2 = a»v 2 / 3 2 + ( 1 - O') n 2 /3 2 

(40) 

m 3 = av 3 /% + ( 1 - a) n 3 /3 3 . 

(41) 

As before, solving equation (39) for v gives 


mi - (1 - a) n/3i 

(42) 

V " OtPi 
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Inserting this expression for v into equation ( 40) and subsequently solving 
for n, we have 


m[(i-a)P 2 (i-a)02 oi[razP\- (m[) 2 /3 2 ] 

n = — ' — — 

( 1-GO01& 

Substituting this expression for n back into (42) gives 

mJ(l-Q!)02 ± n/ a(l-ot)a[m 2 p\ - (mi) 2 /3 2 ] 
a0102 



( 43 ) 


(44) 


Upon substituting the expression for n from (43) and the expression for v 
from (44) into equation (41) , we obtain an equation in the two unknowns 
a and y which may be solved by the iterative process described in the 
general case. We may now solve equations (43) and (44) for 0 2 an< ! 9 j. 

As in the other cases, the quadratic solution in equation (43) results in 
more than one set of estimates. Again, we use the set of estimates which 
gives the closest agreement to the fifth noncentral moment of the sample. 

As before, an alternate solution would be to estimate a graphically 
and then solve the resulting equation (41) directly for y. 

In the event that a is known in advance of sampling we may solve 
equation (41) directly for y and subsequently obtain d 2 and 9 1 from (43) 
and ( 44) . 


Yi = -y 9 = Known . If we let 7i = y% = y, the first three equations of (7) 

become 


mi = adl /y T| 

\y ' 

) + <* 

- a) 

e 2 /y r | 

(H 

(45) 

m 2 = a6^ y T ( 

'* + i' 

<7 i 

) + d 

- a) 

ef /y r | 


(46) 
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(47) 


m 3 


= o,e\ /y 




+ (1 - Oi) 


»? /y 



Thus, we must estimate a, and d 2 only. For simplification, we will ’et 



Now, equations (45), (46) and (47) become 


adl^ + (i - a) = C! 

<xb\^ + ( 1 - a) = c 3 • 

Solving equation ( 48) for a we have 


(48) 

(49) 

(50) 


a = 


c i ~ 0» 


1 / y 


o\ h - o, 171. 


(51) 


Substituting the expression for a from (51) into (49) and (50) , we have after 
considerable algebraic manipulation the equations 


o\ h el /y 


= C, (eY y + «l /y ) - o 2 


(52) 


e} /y o} /y ( 4 /T + » 2 1/T ) = C, (e* /y * o\ ,y e 2 1/y ♦ ef^) - c 2 . (53) 
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1 Ay i Ay 

Inserting the expression for 6 X 0 2 from (52) into the left side of equation 
( 53) , we obtain after simplifying 


J/y 0 Vy 

“l #2 C 1 ” c 2 


( e ;/r + e i/r) +0j 


= 0 


(54) 


Solving (54) for Q^/ 7 gives 


l/ y 

i/y = 0 2 ' c 2 - c 3 
1 .l/y 

$2 Cj - C 2 


(55) 


i/v 

Substituting this expression for 0/ back into equation (52) , we have after 
simplification the quadratic equation in 6. 7 


(c 2 - c\) [e^ 7 ] 2 - (c 3 - CiC 2 ) 6^ y + (c x c 3 - c|) = 0 , (56) 


whose solution for 0^ is 


0 i = 


2(c 2 - c 2 ) 1 £(-c 3 + Cjc 2 ) ± (c 3 - 6cic 2 c 3 - 3c 2 c 2 + 4c 2 c 3 + 4c|) 1//2 


-.r 


• (57) 


Without loss of generality, we may impose the restriction that 0 t < 0 2 . Thus, 
we obtain 9 X and 0 2 from equation (57) using the negative and positive roots, 
respectively. Once we have determined 0! and 0 2 , we obtain our estimate 
for a from equation ( 51) . 

= §2 ~ Known . If we let 6 X = 0 2 = 0, the first three equations of (7) 

become 
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mj = G!0^ yi /3j + (1 - a) 0^ Yz 

(58) 

m 2 = oi0 2/fyi fa + (i - a) 0 2//yz ip 2 

(59) 

= aG 3//ri p 3 + (1 - a) G 3 ^ 1 ip 3 . 

(60) 

Thus , we must estimate the parameters a , yj and y 2 only, 
solving ( 58) for a gives us 

As previously, 

mi - e 1/Y2 i Pi 
“ " 8 1/y ‘ S, - e 1/T * * 

(61) 

Substituting this expression for a into equation ( 59) , we have after simplifi- 
fying 

ml (<, 2/y > h - » 2/y * *,) + mi (e l/ * *, - e 1/y > *,) 


- s 2/y > ft e 1/y ’ *, + e 2/y ’ fc e 1/Tl p, = o . 

(62) 

Equation (62) is an equation in the two unknowns yj and y 2 which may be 
solved by the iterative procedure described for the general case of estimation. 
Once y t and y 2 are determined, we obtain our estimate for a from equation 
(61). 

0i = Go = Unknown. For this special case, it is suggested that we solve 
for a graphically and subsequently follow the procedure outlined for the 
general case of estimation. 

1 
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An illustrative Example 


To illustrate the estimation procedure outlined in this paper for the 
general case, we will consider a sample of 2000 observations selected from a 
mixed population constructed by combining two Weibull distributions with 
■yj = 2. 0000, = 10. 0000, y 2 = 0. 8000, 0 2 = 1. 0000 and a = 0. 8000. The 

sample is summarized in Table XII. For the sample selected, mj = 2. 4708, 
m^ - 8. 6270, m 2 = 36. 3408, mj = 174. 9190 and mj = 935. 3733. 


TABLE XII. A SAMPLE OF 2000 OBSERVATIONS FROM 
A MIXED WEIBULL POPULATION 


CLASSES 

CLASS MARKS 

fl 

h 

f 

CUMULATIVE FREQUENCY 
IN PERCENT 

0 - 0. 5 

0. 25 

40 

175 

215 

10. 75 

0. 5 - 1. 0 

0.75 

113 

78 

191 

20. 30 

1. 0 - 1. 5 

1. 25 

170 

47 

217 

31. 15 

1. 5 - 2. 0 

1. 75 

205 

30 

235 

42. 90 

2. 0 - 2. 5 

2. 25 

216 

20 

236 

54. 70 

2. 5 - 3. 0 

2. 75 

206 

14 

220 

65. 70 

3. 0 - 3. 5 

3. 25 

180 

10 

190 

75. 20 

3. 5 - 4. 0 

3.75 

147 

7 

154 

82. 90 

4. 0 - 4. 5 

4.25 

112 

5 

117 

88. 75 

4. 5 - 5. 0 

4.75 

80 

4 

84 

92. 95 

5. 0 - 5. 5 

5. 25 

54 

3 

57 

95. 80 

5. 5 - 6. 0 

5.75 

34 

2 

36 

97. 60 

6. 0 - 6. 5 

6. 25 

20 

1 

21 

98. 65 

6. 5 - 7. 0 

6.75 

11 

1 

12 

99. 25 

7. 0 - 7. 5 

7.25 

6 

1 

7 

99. 60 

7. 5 - 8. 0 

7.75 

3 

1 

4 

99. 80 

8. 0 - 8. 5 

8.25 

1 

1 

2 

99. 90 

8. 5 - 9. 0 

8.75 

1 


1 

99. 95 

9. 0 - 9. 5 

9.25 

1 


1 

100. 00 



1600 

400 

2000 

_ J 


In Table XU, fj = class frequencies from f^x) , f 2 = class frequencies from 
f 2 (x) , and f = class frequencies from the resulting mixed distribution. 



Figure 4 is a graph of the compound density function and its compo- 
nent distributions. Notice at this point that y Si produces a J-shaped function 
while y > 1 produces a bell- shaped curve. 



FIGURE 4. ILLUSTRATIVE EXAMPLE OF MIXED 
WEIBULL DENSITY FUNCTION 


Employing the graphical technique described in the first section pro- 
vides an estimate of a, the proportionality parameter, equal to 0. 80 as shown 
on Figure 5. 

Once our estimate of a has been determined graphically, we solve 
equation (11) iteratively for first approximations to yi and y 2 . Corresponding 
first approximations for and 0 2 are obtained from equations (13) and (12) , 
respectively. Each set of first approximations is introduced into the fourth 
equation of ( 6) to approximate the fourth noncentral theoretical moment, n\. 
Each set of first approximations is also substituted into the final equation of 
(6) to approximate the fifth noncentral theoretical moment as suggested 
by Pearson [4] . The set of estimates which gives the closest agreement 
between the fifth noncentral moment of the sample mj and the corresponding 
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PERCENT 




"fitted" compound curve given by the final equation of ( 6 ) is the required set of 
estimates. Utilizing the computer program given in the Appendix, we find 
that our estimate of y* lies between 1. 90 and 2. 10. The corresponding 
estimates for the remaining parameters, and M 5 , are as follows: 


7l 

72 


02 

Mi 

Mi 

1. 90 

2. 30 

11. 52 

1. 29 

343. 8875 

1926. 6399 

2 . 10 

0 . 80 

6.74 

0 . 82 

64. 6006 

287. 5141 


Our next approximation to y* is obtained by simple linear interpolation as 
indicated below. 


78 



7i 


1. 90 

343. 8875 

2. 0210 

174. 9190 

2. 10 

64. 6006 


Substituting y t = 2. 0210 into equation (11) , solving for y 2 and subsequently 
solving equations (13) and ( 12) for and 0 2 > we obtain y 2 =0. 8339, 

9 1 = 10. 1496 and 0 2 = 0. 9929. Introducing this set of approximations into 
the final equation of ( 6) gives n\ = 935. 3646. This value of n\ is in such 
close agreement with the corresponding sample moment mg = 935. 3733 that 
we are justified in accepting this set of approximations as our final esti- 
mates. However, if further preciseness is desired, this iterative process 
may be continued to any desired degree of accuracy. 

The computer program outlined in the Appendix gives all possible 
solutions to the estimating equations. Specifically, equations (12) and (13) 
produce four solutions resulting from the four combinations of the positive 
and negative signs prefixing the radicals; i. e. , the combinations are (-,-), 
.(-,+) , (+, -) and (+,+) . It was discovered that the computer program 
printout gave closest agreement between the theoretical moment and the 
sample moment mj when the combination (-,-) was used. 

Comparisons of expected with observed frequencies, along with a com- 
parison of observed and expected distribution functions, are presented in 
Table XIII. It also seems appropriate to compare the observed frequencies 
for the mixed sample with these same frequencies assuming that the sample 
fits a simple Weibull distribution. This will prove or disprove that our mixed 
sample could be treated as a simple Weibull distribution. Cohen's maximum 
likelihood estimation procedure [16], was used to derive estimates for the 
parameters, and the resulting expected frequencies were obtained. Notations 
used in Table XIII are as follows: 

f = observed frequencies for mixed data. 

F = observed distribution function for mixed data, 
o 

f = expected frequencies using estimates derived in this report. 

© 
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F t = expected distribution function using estimates derived in this 
report. 

f = expected frequencies assuming data fits a simple Weibull 
1 3 distribution. 

F = expected distribution function assuming data fits a simple 
eS Weibull distribution. 


TABLE XHI. OBSERVED AND EXPECTED FREQUENCIES FOR 2000 
OBSERVATIONS FROM A MIXED WEIBULL DISTRIBUTION 



F 

es 

0 . 

0717 

0 . 

1932 

0 . 

3292 

0 . 

4621 

0 . 

5821 

0 . 

6843 

0 . 

7677 

0 . 

8332 

0 . 

8828 

0 . 

9195 

0 . 

9458 

0 . 

9642 

0 . 

9768 

0 . 

9852 

0 . 

9908 

0 . 

9943 

0 . 

9966 

0 . 

9980 

0 . 

9988 


The agreement between observed frequencies for the sample and ex- 
pected frequencies using the derived estimates is very good as shown in Table 
XIII. The corresponding observed and expected distribution functions are in 
very close agreement with the maximum absolute difference of 0. 0050 occurring 
at the class 3. 5 - 4. 0. Comparing this value with the Kolmo go rov-Sm i r no v 
statistic , we see that 
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D = D* 0 ®? = 0. 0364 
a o. oi 

gives an. excellent n goo'dness-of-fit" at the 99-percent level of confidence. 

Comparing F with F shows a maximum absolute difference in the 
o es 

distribution functions of 0. 0358 occurring at the class 0 - 0. 5. Since the 
Kolmogorov-Smirnov statistic 

D N = = 0. 0304 , 

CL 0-05 

this value of 0. 0358 is sufficient to reject the hypothesis that our mixed distri- 
bution could be considered a simple Weibull distribution. 

As an alternate goodness-of-fit test for agreement between observed 
frequencies and expected frequencies using the derived estimates, the x 2 
index was calculated and the results are as follows: 

X 2 = 1. 2669 d. f. = 10 P(X 2 > 1. 2669 = 0. 995) . 

Thus, in consideration of the low x 2 index of dispersion,, we may conclude 
that we have an excellent fit for the chosen sample; 

Conclusions 

It is an accepted fact that the method of moments is not (except for 
distributions such as the normal, binomial, and Poisson) the most efficient 
procedure for estimating the parameters of a frequency distribution. Methods 
having maximum efficiency, such as the method of maximum likelihood, are 
more desirable. However,- in the case of the mixed Weibull distribution with 
its five parameters, the maximum likelihood estimating equations are almost 
intractable. 

The central, noncentral, and factorial moments of this distribution 
were investigated, and it was discovered that the noncentral moments possessed 
optimum characteristics for the development of estimating equations. A 
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comparison of noncentral sample moments with the theoretical moments for the 
sample selected showed an error of 0. 09 percent for the first moment, 0. 44 
percent for the second moment, 1. 64 percent for the third moment, 4. 94 
percent for the fourth moment, and 12. 73 percent for the fifth moment. This 
progressive increasing percentage of nonagreement between sample moments 
and theoretical moments illustrates the large sampling errors involved in the 
use of higher order moments. Using only the first three sample moments 
with their relatively low percentage of error in the estimating equations pro- 
duced very good agreement between final estimates and the population 
parameters yj, y 2 , 0* and 0 2 . Sheppard's corrections for grouped data were 
applied to the sample moments in order to increase this agreement, but pro- 
duced no significant change in the results; therefore, the corrections were not 
used in the estimating equations. 

This paper presents an estimating procedure that produced very good 
results for the sample chosen. The use of electronic digital computers makes 
the somewhat involved method practical and applicable to experiments in which 
the mixed Weibull distribution is the appropriate statistical model. 

In the author's opinion, the estimating procedures outlined in this 
report warrant further investigation for increase in efficiency, and improve- 
ment and simplification of form for the estimating equations involved. 


George C, Marshall Space Flight Center 

National Aeronautics and Space Administration 
Huntsville, Alabama, March 29, 1967. 
160-44-04-00-62 
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Appendix 

Computer Program for Estimating the Parameters 
of a Mixed Weibull Distribution 


3200 FORTRAN 

ESTIMATION IN MIXTURES OF TWO WEIBULL DISTRIBUTIONS 


n ! MFNS T ON XT 1*0) »F I f 30 ) 

RFAOl A0*1 )DG1 *0G?*FG1 *PG? q) 

RFAD(60»1 )PGAM1 *RGAM?,FN*ALP 
1 FORMAT ( 4F1 0 . * ) 

RPAH( AO*? )NV1 *K 
7 FORMAT ( 13/13) 

RFADt AO** ) *X| ( T ) *T = 1 *NVl 1 
R F AD ( AO*** J (FT IT 1 * !*1 *NV1 1 
3 FORMAT ( A^l 0*e ) 


VIP I TP ( A1 ,A ) 

13 

PM 1 =0.0 


PM?=0.0 


DM1=0,0 

14 

PM4=0 . 0 


pM A =0 . 0 


on «; t = i,mvi 

13 

pr 1 =X T ( T > 


FR? = FP7 #XJ ( T ) 


PR 3 = FR ?* X I ( I ) 

16 

PR4=PR3*XT 1 T ) 


PR S = FR 4*X I ! I ) 


PM1 =PM1+FR1*FT ( I ) 

1 7 

PM? = PM? + PP2*FT < I ) 

27 

PM3 = PM3+PR3*F T ( I ) 


PM4=PM4+pP4*FT ( T 1 

7 ft 

omf = PM^ + pd«;*f| ( 1 > 

79 

rnNT! NU p 


PM1=PM1/FM 


p M?=PM2/FN 

1 0 


PM3=PM3/FN 

PM4=PMA/* r N 

pMA=PMA/PM 

WPI TP ( A1 *70 »PM1 ,pM?,PM1 ,OM4 ,PMP 

G AM 1 =PGAM1 

GAM2=PGAM? 

23 TRM=1 ,0/GAMl+l .0 
ni=GAMMA(TRM) 

TRM=?.0/GAM1+1 .0 
R?=GAMMA(TRM) 

TRM=3,0/r-AMl + 1 .0 
R^aAAMMA { TRM ) 

WR f TP ( A 1 *o?) 

O0 TRM=1 .0/GAM7+1 ,0 
PS I 1 = G AM M A < TRM) 

TRM=?.0/GAM?+1 .0 
PSI?=GAMMAtTPM) 

TRM=3,0/GAM?+1 # 0 
OSI3=GAMma ( TRM) 

FR = ALP* ( 1 .0-ALP)*Rl #*7#PS I 7 
PR 1 = ( 1 .0-ALP )**?*PSI1**?*R? 
DFN=?.0*(FR+FP1 ) 
pR7=-4,0*PM1 ** 7 *FR *R ? 

FR3=A.0*PM2*ALP»FR1*B1 **? 

PRA = 4.0*PM?»AIP**?*( 1.0-AI.P)*P1 **4*PSt? 
TRM=PR?+PR3+FP4 


r F ( TPM) 17*91 ♦ R 1 
TPM1=S0RT(TPM) 

FRA=?.0*PM1 *( 1 .0-ALP)*PSIl*B? 

FNP=( FRA+TRM1 )/DFN 
PNM=T PR5-TRM1 )/OFN 
ER6=( PMl-d ,0-ALP)*PSll ) /(ALP*R1 ) 

VP=FRA*FNP 
VM=ER6*FNM 
GO TO ( 1 7,1 4*1 3,1 A>* 

V = VP 
FN=FNP 
GO TO 10 

V = VP 
fn=fnm 
GO TO 10 

V = VM 
FN=ENP 
GO TO 10 

V = VM 
EN=FNM 
GO TO 10 

1 P{ FG?-OAM?-.000001 )?R.?R»?7 
GAM2=GAM?+DG? 

GO TO 90 

IF(FG1-GAM1_.P0001 ) 1 90, 1 oo,?q 
GAM1 =GAM1 +HG1 
GAM2=RGAM? 

GO TO 23 

TH1 =FXPCGAM1*AL0G( V) 1 
TH2=FXP(GAM2*ALOG(FN) ) 

ER=S,0/GAM1+1.0 

TRV= GAMMA ( FR ) 

PR = A # 0/GAM24.1 .0 
TRMi=r,AMMA(FP> 
pp 1 = AL p* ( TH1 **< 3.0/GAM) ) ) 

FR?= ( 1 .O-ALP ) *( TH?**( S. O/GAM? ) ) 
UP5=FR1*TRM+FR2*TRM1 
WRITF(61*36»GAM1 ,GAM2,TH1 ,TH2*UPS 
GO TO 12 
100 STOP 

4 FORMAT ( 1H1 ) 

70 FORMAT {3X,AHP M 1 =Fl 4 . 7 , 3X , AHPM? = F 1 4 .7 , 3X , AHOW|3 
1 3 X , AHPM4 =F1 4*7,7X ,AHPM«> =F14.7/> 

92 FORMAT!//) 

26 F0RMAT(3X,6HGAM1 =F1 3 . * , *X , 6HGAM2 =F13.A,3X, 

16HTH1 =F13.S ,3X,6HTH2 =E 1 3 . P , 3X , AHU5 =E13.«) 

END 


=r l 4. 7 , 
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^UNCTION GAMMA (X) 

7R FORMAT ( RRH GAMMA FUNCTION OF A MFGATTVF INTFGFR IS NOT OFFIN^D.) 

5 I F ( X ) 10 ♦ 24 ♦ 1 

10 N = -X 
EN=-N-1 
V=X-EN 

I F C V— 1 * ) 70 » 80, 2 0 
1 *4 N = X 
EN=N 
V=X-FM 

70 GAMMA=1 ,+V*( .A??7BA^A +V*(.411»407R + V* ( .OR 1 *7871 9 + V* 

l(.074?77o08 +V*(-. 00071 0007A7 +V* (. 01007? 696 + V* 0074 6674 R 

?+V*( .001 RT076fli -V*( ,000'»A47RA70 -V* . 0000677 1 0*7 ))))))))) 

I F ( FN-? • ) 37 » 25 * BO 

24 'GAMMA* 1« 

25 RETURN 
30 N=N-.l 

DO 35 I = 2 * N 
F 1 = I 

35 GAMMA=GAMMA*(FI+V) 

RFTURN 
37 N=2.-FN 

DO 40 1=1 »N 
FI =2-1 

40 GAMMA=GAmma/(F1+V) 

RFTURN 

BO WRITE(61»75) 

STOP 

END 


GAM00001 

GAM00007 

GAM00003 

GAM00004 

GAM00005 

71M00006 

GAM00007 

GAM00008 

GAM00009 

GAM00010 

GAM00011 

GAM0001 7 

GAM00013 

GAM00014 

GAM00015 

GAM00016 

GAM00017 

GAM00018 

GAM00019 

GAM00020 

GAM00021 

GAM00022 

GAM00023 

GAM00024 

GAM00025 

GAM00026 

C5107570 

GAM00030 


FINIS 



NOTE: TRM IS THE QUANTITY UNDER THE RADICAL IN EQUATIONS 12 AND 13. 

VP IS THE VALUE OF EQUATION 13 WHEN THE POSITIVE VALUE IS USED. 
VM IS THE VALUE OF EQUATION 13 WHEN THE NEGATIVE VALUE IS USED. 
ENP IS THE VALUE OF EQUATION 12 WHEN THE POSITIVE VALUE IS USED. 
ENM IS THE VALUE OF EQUATION 12 WHEN THE NEGATIVE VALUE IS USED. 


COMPUTER PROGRAM FLOW CHART FOR ESTIMATING THE PARAMETERS OF A MIXED WEIBULL DISTRIBUTION 


00 

01 






Load Data 


First, two tables are read in starting in XI and FI 

DG1 - Ayi 

DG2 - Ay 2 

EG1 - Terminal y t 

EG2 - Terminal y 2 

BGAMi - Beginning yi 

BGAM2 - Beginning y 2 

EN - Sample size. 2000 used for this program 

ALP - a 0. 8 used 

NVI - Number of values in the tables 

k - An integer 1 4 

The computed Ms should approach mj 

This can be accomplished by a series of runs varying Ay t , Ay 2 , y^ and y 2 
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Legend 


mi - PMi e t - THl 

mi - PM 2 0 2 - TH2 

mj - PM 3 
mi - PM 4 
mi - PM5 

01 ~ Bi 

02 - B2 

03 " B3 

! pi - PSIi 
ip 2 - PSI2 
ip 3 - PSI3 
Ms -UP5 

N - EN Sample Size 

n - EN Once the above N is used the location EN is no longer needed. 

It is used to store n. 

v - V 

ENP - In equation ( 12) , if the positive value of the radical is used. 

ENM - If the negative value is used. 

VP - In equation ( 13) , if the positive value of the radical is used. 

VM - If the negative value is used. 

The following four combinations are possible, depending upon what number is 
read into location k. 

v n k Regardless of which combination is used, before the 

+ + 1 program is continued, the value of equation (13) is 

+ - 2 stored in V and the value of equation ( 12) is stored 

+ 3 in EN. 

4 
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